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There have been many attempts to derive continuum models for dense granular flow, but a 
general theory is still lacking. Here, we start with Mohr-Coulomb plasticity for quasi-2D granular 
materials to calculate (average) stresses and slip planes, but we propose a "stochastic flow rule" 
(SFR) to replace the principle of coaxiality in classical plasticity. The SFR takes into account two 
crucial features of granular materials - discreteness and randomness - via diffusing "spots" of local 
fluidization, which act as carriers of plasticity. We postulate that spots perform random walks biased 
along slip-lines with a drift direction determined by the stress imbalance upon a local switch from 
static to dynamic friction. In the continuum limit (based on a Fokker-Planck equation for the spot 
concentration) , this simple model is able to predict a variety of granular flow profiles in flat-bottom 
silos, annular Couette cells, flowing heaps, and plate- dragging experiments - with essentially no 
fitting parameters - although it is only expected to function where material is at incipient failure 
and slip-lines are inadmissible. For special cases of admissible slip-lines, such as plate dragging 
under a heavy load or flow down an inclined plane, we postulate a transition to rate-dependent 
Bagnold rheology, where flow occurs by sliding shear planes. With different yield criteria, the SFR 
provides a general framework for multiscale modeling of plasticity in amorphous materials, cycling 
between continuum limit-state stress calculations, meso-scale spot random walks, and microscopic 
particle relaxation. 



I. INTRODUCTION 

For centuries, engineers have described granular ma- 
terials using continuum solid mechanics [J Q • Dense 
granular materials behave like rigid solids at rest, and 
yet are easily set into liquid-like, quasi- steady motion 
by gravity or moving boundaries, so the classical the- 
ory is Mohr-Coulomb plasticity (MCP), which assumes 
a frictional yield criterion. The simplest model is the 
two-dimensional "Ideal Coulomb Material" at limit-state, 
where the maximum ratio of shear to normal stress is ev- 
erywhere equal to a constant (the internal friction coef- 
ficient), whether or not flow is occurring. This model is 
believed to describe stresses well in static or flowing gran- 
ular materials, but, as we explain below, it fails to pre- 
dict flow profiles, when combined with the usual Coaxial 
Flow Rule of continuum plasticity. Indeed, it seems con- 
tinuum mechanics has not yet produced a simple and 
robust model for granular flow. 

In recent years, the sense that there is new physics 
to be discovered has attracted a growing community of 
physicists to the study of granular materials 4, 5, |y, 0, 
Unlike the engineers, their interest is mostly at the 
discrete particle level, motivated by the breakdown of 
classical statistical mechanics and hydrodynamics due to 
strong dissipation and long-lasting, frictional contact net- 
works. Dense granular materials exhibit many interesting 
collective phenomena, such as force chains, slow struc- 
tural relaxation, and jamming. Similar non-equilibrium 
phenomena occur in glasses, foams, and emulsions, as 
in granular materials, so it is hoped that a general new 
statistical theory may emerge. Presumably from such 
a microscopic basis, continuum models of glassy relax- 
ation and dense granular flow could be systematically 
derived, just as dissipative hydrodynamics for granular 



gases can be derived from kinetic theory with inelastic 
collisions 0. 

This dream has not yet been achieved, but many em- 
pirical continuum models have been proposed 0, l^, El • 
The difficulty in describing dense granular flow is evi- 
denced by the remarkable diversity of physical postulates, 
which include: coupled static and rolling phases 

m Hi m 

jBagnold rheology based on "granu- 
lar eddies" [T3, granular temperature-dependent viscos- 
ity [3 5 density-dependent viscosity 20], non-local 
stress propagation along arches 0, self-activated shear 
events due to non-local stress fluctuations [2J. [2J. free- 
volume diffusion opposing gravity 24, Ig^ l26ll?n. [2^ . 
"shear transformation zones" coupled to free- volume ki- 
netics |2^ I33, and partial fluidization governed by a 
Landau-like order parameter |^ [s^ • Each of these the- 
ories can fit a subset of the experimental data 33], usu- 
ally only for a specific geometry for which it was de- 
signed, such as a flowing surface layer 11^ [T^ 113. ITH l3lj . 
inclined plane 0, Couette cell *19|, l20l . inclined 
chute 22, 23], or wide silo 24, 25, 2&,^l2j and none 
seems to have very broad applicability. For example, we 
are not aware of a single model, from physics or engineer- 
ing, which can predict velocity profiles in both draining 
silos and annular Couette cells, even qualitatively. 

The theory of partial fluidization of Aranson and Tsim- 
ring has arguably had the most success in describing mul- 
tiple flows within a single theoretical framework 0, • 
Although setting boundary conditions for the order pa- 
rameter usually requires additional ad hoc assertions, the 
model is nonetheless able to reproduce known flow behav- 
ior in inclined chutes, avalanches, rotating drums, and 
simple shear cells without many fitting parameters. It 
also describes some unsteady flows. However, the theory 
lacks any clear microscopic foundation and is not directly 
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coupled to a constitutive stress model for static materials. 
As such, it has only been applied to problems with very 
simple solid stress fields, limiting its current applicability 
to fiows that depend on only one spatial variable. 

In an attempt to describe arbitrary geometries, such as 
silos and Couette cells, we take the view that the engi- 
neers may already have a reasonable continuum descrip- 
tion of the mean stresses, so we start with Mohr-Coulomb 
plasticity. However, discreteness and randomness clearly 
need to be taken into account in a granular material. For 
static stresses, quenched randomness in material proper- 
ties is known to lead to statistical slip-line blurring in 
"stochastic plasticity" 34], but this says nothing about 
how plastic yielding actually occurs. 

To describe yielding dynamics^ we propose a "stochas- 
tic fiow rule" (SFR) where local fiuidization (stick-slip 
transition) propagates randomly along blurred slip-lines. 
We build on the recently proposed Spot Model for 
random-packing dynamics 27] by viewing "spots" of free 
volume as carriers of plasticity in granular materials, 
analogous to dislocations in crystals. Multiscale spot 
simulations can reproduce quite realistic fiowing packings 
in silo drainage |2^; here, we introduce a mechanical ba- 
sis for spot motion from MCP, which leads to a theory 
of considerable generality for bulk granular fiows. 

The paper is organized as follows. Since plasticity is 
unfamiliar to most physicists, we begin by reviewing key 
concepts from MCP in section Ull both for stresses and 
for dense fiows. In section Hill we highlight various short- 
comings of the classical theory, many of which we at- 
tribute to the Coaxial Flow Rule. We then introduce the 
general spot-based SFR and a specific simplification to 
be used for granular fiow in section HVl Next we apply 
the theory to four prototypical examples: silo, Couette, 
heap, and plate-dragging fiows in section Then in sec- 
tion we explain how the last two examples indicate 
a smooth transition from the SFR to Bagnold rheology, 
when slip-lines become admissible, and we present a sim- 
ple composite theory, which extends the applicability of 
the model to various shear fiows. In section IVTll we con- 
clude by further clarifying the range of applicability of 
the SFR and possible extensions to other granular fiows 
and different materials. 



II. CONCEPTS FROM CONTINUUM 
MECHANICS 

A. Mohr-Coulomb plasticity: stresses 

In the eighteenth century, it was Coulomb, as a mil- 
itary engineer designing earthen fortresses, who intro- 
duced the classical model of a granular material, which 
persists to the present day: a continuous medium with 
a frictional yield criterion. His ideas were expressed in 
general continuum-mechanical terms by Mohr a century 
later, and a modern mathematical formulation of "Mohr- 
Coulomb plasticity" (MCP), which we also use below, is 




FIG. 1: Stresses on a material element. All vectors are point- 
ing in the positive direction as per our sign convention. 
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FIG. 2: Force diagram for a wedge. Hypotenuse length as- 
signed to unity. 



due to Sokolovskii0. Although other mechanical mod- 
els exist, such as Drucker-Prager plasticity |^, MCP is 
perhaps the simplest and most widely used for granular 
materials in engineering 1]. As such, we choose to build 
our model of dense granular fiow on the MCP description 
of stresses, as a reasonable and time- tested first approx- 
imation. 

We begin in this section by reviewing relevant con- 
cepts from MCP, e.g. following Nedderman 1]. The 
fundamental assumption is that a granular material can 
be treated as an "Ideal Coulomb Material" (ICM), i.e. a 
rigid-plastic continuous media which yields according to 
a Coulomb yield criterion 



\r/(j\ =11 = tanc 



(1) 



where r is the shear stress, a is normal stress, and (j) 
the internal friction angle, akin to a standard friction 
law with no cohesion. Throughout, we accept the com- 
mon tensorial conventions for stresses with the key ex- 
ception that normal stresses are deemed positive in com- 
pression. This is a standard modification in the study 
of non-cohesive granular materials since granular assem- 
blies cannot support tension. We will also focus entirely 
on quasi-2D geometries. 

Consider a small material element in static equilibrium 
and with no body forces present (see Figure QJ. The nor- 
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mal stresses a. 



and (J yy can differ and the shear stresses 



xy and Tyx must be equal in order to balance moments. 
Likewise the variable Tyx is redundant and will not be 
used again in this paper. To determine the stresses along 
any angle within this element, we place a new boundary 
within the material at some desired angle and observe 
force balance on the wedge that remains (see Figure [2)). 
After algebraic simplification, this gives 

^6* = 2 (^^^ + ^2/2/ ) + 2 '^^^'^ ~ ^ ~ ^""y 

= ^ {axx - cTyy ) sin 26 + Txy cos 26 



Now define 



tan 2^|J 
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R = 



which allows us to write 



''XX ^ yy 



ae = p + Rcos{2e -2ip) 
Te = Rsm{2e-2iP) 



(2) 
(3) 



This implies that for all angles 6>, the locus of traction 
stresses {ere ^ re) is a circle centered at (p, 0) with radius 
R. This circle is referred to as "Mohr's Circle" . 

We have just derived Mohr's Circle without accounting 
for the possible effects of body forces acting on the mate- 
rial element and gradients in the stress field. Adjusting 
for these effects, however, would change the results only 
negligibly as the element gets small in size. If we were to 
apply the same force-balancing analysis to a differentially 
small material element with a body force and stress gra- 
dients, we would find that the stress differences on the 
walls and the inclusion of the differentially small body 
force within only add differentially small terms to the 
equations for ae and re. Thus we can always use Mohr's 
Circle to obtain traction stresses along a desired angle. 

To ultimately define a full stress state for the material 
element, we need one more equation — we have 3 stress 
variables and only 2 force balance equations: 



(4) 
(5) 



dcTxx 

dx 


drxy 
dy 


rpx 

body 


ddyy 

dy 


drxy 
dx 


- py 

body 



We say a material element is at incipient failure if 
the yield criterion is fulfilled along some direction and 
|T/cr| < /i along all others. A material in which incipient 
failure occurs everwhere is said to be at a limit-state. In 
a limit-state, the Mohr's Circle at every point in the ma- 
terial must be tangent to the locus \r/a\ = fi. As can be 



seen by applying trigonometry in Figure Ol this require- 
ment means that R = p sin (/), enabling us to parameter- 
ize the stresses in terms of p and t/j only, thereby closing 
the equations. For this reason, we restrain our anlysis to 
limit-state materials and refer to p and as the stress pa- 
rameters or Sokolovskii variables. (The limit-state stress 
treatment described here is also known as "Slip-Line The- 
ory"; to avoid possible confusion, we specify this is not 
equivalent to Limit Analysis Plasticity concerned with 
upper and lower collapse limits.) 

Solving for the original stress variables in terms of the 
stress parameters gives: 



O'x 



yy 



'xy 



= p(l + sin^cos 2?/^) 
= p{l — sin (j) cos 2?/^) 
= —p sin (j) sin 2?/^ 



(6) 
(7) 
(8) 



Using these expressions, we re- write equations (0} and 
©: 



sin (j) sin 2?/^ py 
- sin (j) cos 2tlj)py 



(l+sin0cos2?/^)pa; — 2p sin ^ sin 2?/^ tpx ^ 

+ 2p sin cos 2iIj i/jy = F^^^y 
sincj) sin 2ip px + 2p sin (j) cos 2tlj t/jx -\- {1 
+ 2p sin (/) sin 2?/^ ijjy = py^^y 

These will be referred to as the "stress balance equa- 
tions". They form a hyperbolic system and thus can be 
solved using the method of characteristics. The system 
reduces to the following two characteristic equations: 

dp ^ 2pii dip = py^^y {dy ^ ji dx) ^ F^^^y{dx ± ji dy) 

dy 

along curves fulfilling -— = tan(?/^ ^ e). (9) 
dx 

To solve the stress balance equations, mesh the two fam- 
ilies of characteristic curves in the bulk, then march from 
the boundaries in, progressively applying the two differ- 
ential relationships above to approximate the stress pa- 
rameters at each intersection point in the mesh. More on 
this can be found in 36]. Other ways to solve the stress 
balance equations include the Two- Step Lax- Wendroff 
Method 37] and the Galerkin Method "H^. 

We return now to Mohr's Circle for a discussion of the 
stress properties within a differential material element. 
Equations Q and show that Mohr's Circle can be 
used as a slide-rule to determine the stresses along any 
angle 6: One arrives at the point {ere, re) by starting 
at {(rxxi^xy) and traveling anti-clockwise around Mohr's 
Circle for 26 radians (see Figure Also note on the 
diagram that the stresses along the x and y directions 
lie along a diameter of Mohr's Circle; any two material 
directions differing by an angle of 7r/2 lie along a diameter 
of the corresponding Mohr's Circle diagram. Utilizing 
this property in reverse is perhaps the easiest way to 
draw Mohr's Circle in the first place; draw the unique 
circle for which {(Jxxi^xy) and {(Jyy, ~'^xy) are endpoints 
of a diameter. 

Let (cri,0) and (cr3,0) be the points of intersection 
between Mohr's Circle and the cr-axis, where cti > a^. 



FIG. 3: Using Mohr's circle jointly with the Coulomb internal yield locus (r — ±iJ.a) to determine the traction stresses along 
any plane within a material element. 



These points correspond to the two lines within a mate- 
rial element along which the shear stress vanishes and the 
normal stress is maximal or minimal, ai {as) is called 
the major (minor) principal stress and the line on which 
it acts is called the major (minor) principal plane. 

Mohr's Circle shows that the major principal plane 
occurs at an angle anti-clockwise from the vertical (see 
Figure 01). Thus the major principal stress points along 
an angle tjj anti-clockwise from the horizontal. This is the 
standard physical interpretation of t/j. One might think 
of ip as the angle from the horizontal along which a force 
chain would be predicted to lie. 

By right-triangle geometry, a line segment connecting 
the center of Mohr's Circle to a point of tangency with 
the internal yield locus would make an angle of 7r/2 — ^ 
with the cr-axis. Each point of tangency represents a 
direction along which the yield criterion is met, i.e. a 
slip-line. Mohr's Circle indicates that the slip-lines are 
angled (7r/2 — 0)/2 up and down from the minor prin- 
cipal plane. But since the major and minor principal 
planes are orthogonal, the major principal stress points 
along the minor principal plane. Defining e = 7r/4 — (/)/2, 
we deduce that slip-lines occur along the angles ± e 
measured anti-clockwise from the from the horizontal. 
Looking back at the characteristic equations, we see that 
the slip-lines and the characteristic curves coincide. This 
means that information from the boundary conditions 
propagates along the slip-lines to form a full solution to 
the stress balance equations. 

It is worth noting that the stress balance equations 
are written for static materials and do not appear to ac- 
count for dynamic behavior like dilatancy and convection 
stresses. The theory of critical state soil mechanics 
was the first to rigorously approach the issue of dilatancy 
(see appendix). It concludes that when material attains 
a flow state in which the density field stops changing in 
time, all points in the flow lie along a critical state line of 
the form \r/(j\ = 5 for 5 constant. Since this exactly mir- 
rors the Coulomb yield criterion, we can keep the stress 




FIG. 4: Important lines intersecting each material point: (a) 
Major principal plane / Minor principal stress direction; (b) 
Minor principal plane / Major principal stress direction; (c)- 
(d) Slip-lines. 

balance equations and utilize (5 = /i (as in ^3). As for 
convection, adding the pu ■ Vu term into the stress equa- 
tions couples the stresses to the velocity and makes the 
problem very difficult to solve. The practice of ignoring 
convection is justified by our slow- flow requirement and 
is commonly used and validated in basic solid mechanics 
literature |^ . So we conclude that dynamic effects 
in flowing materials do not preclude the use of the stress 
balance equations in slow, steady flows. 

B. Mohr-Coulomb plasticity: flow rules 

To calculate flow, we assert incompressibility and a 
flow rule — the flow rule is a constitutive law chosen to 
reflect the general behavior of the material at hand. The 
continuous nature of the ICM assumption suggests that 
symmetry should be kept with respect to the principal 
stress planes. Based on this, Jenike proposed adopt- 
ing the Coaxial Flow Rule. The principle of coaxiality 
claims that material should flow by extending along the 
minor principal stress direction and contracting along the 
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FIG. 5: Sketch of the Coaxial Flow Rule. 



major principal stress direction; the principal planes of 
stress are aligned with the principal planes of strain-rate. 
The intuition for this constitutive rule is shown in Fig- 
ure [5J Mathematically, this means that in a reference 
frame where the minor and major principal stress direc- 
tions are the basis, the strain-rate tensor should have no 
off-diagonal components, i.e. 



R^ER^ is diagonal. 



(10) 



where R^ rotates anti-clockwise by ip and E is the strain- 
rate tensor 



E=i(Vu+(Vur). 



(11) 



where u = {u^v) is the velocity. Calculating the (1,2) 
component of the matrix in equation ([Tn|) and setting it 
to zero gives the equation of coaxiality, 

du ^ dv fdu 2^p (12) 

dy dx \dx dy J 

This flow rule has played a dominant role in the develop- 
ment of continuum plasticity theory and will be closely 
analyzed in this work. 

Coaxiality with incompressibility comprises another 
hyperbolic system of equations enabling the velocity field 
to be solved via characteristics: 



du^ tan(?/^ ^ 7r/4) d'u = 



dy 



along curves fulfilling -— = tan(?/^ ^ 7r/4). (13) 

dx 

So, given ilj{x,y) from the stress balance equations, infor- 
mation about the flow travels from the boundaries into 
the bulk along curves rotated 7r/4 off from the princi- 
pal stress planes — using Mohr's Circle, we observe that 
these are the lines for which the shear stress is maximal 
(and the normal stresses are equal). 

Other flow rules have been suggested instead of coaxi- 
laity. Of specific note, A.J.M. Spencer 41] has proposed 
the double-shearing flow rule. Unlike coaxiality which 
can be understood as a simultaneous equal shearing along 
both slip-line families, double-shearing allows the shear- 
ing motion to be unequally distributed between the two 
families in such a way that the flow remains isochoric. 



For steady flows, the double-shearing flow rule is 

dy J \dx dy J 

2u-V?A) (14) 



sin 2tlj 



smc 



dx 
dv du 
dx dy 



It can be seen that when the material neighboring a 
particle rotates in sync with the rotation of the prin- 
cipal planes (i.e. as tracked by the rate of change of V^), 
the right side goes to zero and the rule matches coaxial- 
ity. Under double-shearing, the characteristics of stress 
and velocity align, easing many aspects of the numerics. 
Some recent implementations of granular plasticity have 
utilized principles of double-shearing 42] . Though in this 
paper we deal primarily with the comparison of coaxial- 
ity to our new theory, this equation will be mentioned 
again in a later section. 



C. The rate-independence concept 

We now more fully address the conceptual basis for the 
flow theory just introduced. The theory is fundamen- 
tally different from traditional fluids where force-balance 
(including convection and viscous stresses in the case of 
Newtonian fluids) can be used alongside incompressibility 
to fully determine the fluid velocity and pressure fields. 
Unlike a fluid, granular materials can support a static 
shear stress and thus force-balance plus incompressibility 
alone is an underconstrained system. Rather, the stress- 
strain relationship for granular material is presumed to 
be rate-independent in the slow, quasi-static regime we 
study. 

This concept is best understood tensorially. We can 
rewrite the equations of coaxiality and incompressibility 
equivalent ly as: 



where 



E = ATo 



Stress tensor : 



^ XX 

^xy 



(15) 



(16) 



To = T — -(trT)I = Deviatoric stress tensor, (17) 

and A is a multiplier which can vary in space. Equa- 
tion ([T^ is merely the ratio of the (1,2) component and 
the difference of the (2,2) and (1,1) components of equa- 
tion (|T^ . thus cancelling A, and incompressibility is au- 
tomatic since we relate to the deviatoric stress tensor. 
Equation (|T^ gives a simple and highly general form for 
plastic material deformation applicable to a broad range 
of deformable materials and so it is ideal for illustrating 
the role of rate-dependency. In MCP, we solve for Tq a 
priori from the stress balance equations. A adds the extra 
degree of freedom necessary to make sure the strain-rate 
field is actually compatible with a real velocity field — 
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A is not any specific function of the stress or strain-rate 
variables and it adjusts to fit different velocity bound- 
ary conditions. Thus the stress alone does not imply the 
strain-rate and vice versa. 

Supposing on the other hand that we were dealing with 
a rate- dependent (i.e. visco-plastic) material like Newto- 
nian fluid, the above tensorial equation would still apply 
but we cannot claim to know Tq in advance since ma- 
terial motion changes the stresses. Instead we prescribe 
a functional form for A, like A = viscosity"^ = constant, 
and write the force balance equations in terms of E. Thus 
E is computed very differently for the two cases: in the 
rate-independent case, (P3|) is solved using a known form 
for To and in the rate-dependent case, (P3|) is solved us- 
ing a known form for A. 

The physical intuition for rate-independent flow can 
be easily understood with an example. Suppose we slide 
two frictional blocks against each other at two different 
non-zero sliding rates. In most rudimentry dry friction 
laws, the shear stress required to slide one block against 
another is proportional to the normal stress — there is no 
mention whatsoever of the rate of sliding. Thus the two 
sliding rates are modeled to be attainable with the same 
shear stress and likewise the stress-strain relationship is 
deemed rate-independent. For slow granular flows with 
long-lasting interparticle contacts, comparisons with this 
example are especially instructive. 



III. SHORTCOMINGS OF MOHR-COULOMB 
PLASTICITY 

The use of the stress balance equations with incom- 
pressibility and the Coaxial Flow Rule will be referred 
hitherto as Mohr-Coulomb Plasticity (MCP). The the- 
ory has the benefit of being founded on mechanical prin- 
ciples, but does have some marked drawbacks. We point 
out a few: 

• The theory frequently predicts highly discontinuous 
velocity fields. 

• The Coaxial Flow Rule is conceptually troubling in 
some simple geometries. 

• The assumption of limit-state stresses is overreach- 
ing. 

• MCP is a continuum theory and thus cannot model 
discreteness and randomness. 

We will now discuss these four points in detail. 



A. Discontinuous "shocks" in stress and velocity 

The two stress PDEs and two flow PDEs are each fully 
hyperbolic systems meaning that continuous solutions do 
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stress Discontinuity 



Side B 



FIG. 6: Stresses on a control volume intersected by a discon- 
tinuity. Note how a jump in a± places no net force on the 
control volume. 



not necessarily exist for arbitrary choices of the bound- 
ary conditions. Instead, discontinuous solutions are con- 
structed utilizing intuitive jump conditions. For stresses, 
a jump in the stress parameters across a discontinuity line 
is only allowable if such a jump places no net forces on a 
small control volume surrounding the line thereby ensur- 
ing particle stability. This means the normal and shear 
stresses along the direction of the discontinuity must be 
the same on both sides of the discontinuity. However, the 
normal stress along the perpendicular direction can have 
a jump upon crossing the discontinuity as this places no 
net force on the control volume (see Figure Ej). 

In terms of the stress parameters, this means that p 
and ip can jump as long as 

(1 + sin0cos(26 - 2iIj^)) _ sin(2e - 2^/^^) _ 



(1 + sin (j) cos(2e - 2V^^)) sin(2e - 2?/^^) 



p- 



(18) 

where O is the angle from the vertical along which the 
stress discontinuity lies. 

As for velocity, incompressibility forces us to impose a 
simpler jump rule in that the component tangent to the 
velocity discontinuity is the only one allowed to jump. 
We note that whenever a stress shock exists, the jump in 
the stress parameters will usually place a jump in the flow 
rule and may cause a velocity discontinuity to form coin- 
cident with the stress shock. A velocity discontinuity can 
form even when the stress field is smooth since the veloc- 
ity PDEs are themselves hyperbolic. (It is worth noting 
that when shocks are allowed in the solution, multiple 
solutions sometimes arise to the same problem; introduc- 
tion of the so-called "entropy condition" can be used to 
choose the best of the possible solutions 0, Over- 
all, the MCP equations are mathematically very poorly 
behaved, and have also been shown to give violent insta- 
bilities and finite-time singularities ^J, |43 • 

Aside from its mathematical difficulties, MCP theory 
also does not match experiments or our everyday expe- 
rience of granular flows. In particular, MCP commonly 
predicts complicated patterns of velocity discontinuities 
in situations where experiments indicate smooth flow in 
steady-state. In Figure [3 the numerically determined 
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FIG. 7: Numerical solution to MCP in a wedge hopper with 
non-radial stresses on the top boundary. Normal stress in 
the radial direction displayed. (Courtesy of the authors of 
Ref. ^\.) 



stress field for a wedge hopper with only slightly non- 
radial boundary conditions on the top surface exhibits 
a fan-like array of shocks. The associated velocity field 
(not shown) will at best exhibit a similar pattern of dis- 
continuities and at worse add even more discontinuities. 
Such a broken velocity field is clearly unphysical. As the 
grain size becomes very small (sands), discontinuous ve- 
locity fields with no relationship to the stress field have 
been observed, but these are only temporary; the shocks 
commonly blur away immediately after the onset of flow, 
which has been attributed to some instability mechanism 
|46| . Literature on the topic [Ij is quick to concede that 
infinitessimally sharp velocity jumps are physically non- 
sensical and should be understood as being spread over at 
least a few particle widths. Below, we will see that our 
model naturally provides a mechanism for the blurring 
of velocity shocks even in the presence of a stress shock, 
with large velocity variations occurring only at the scale 
of several particle diameters. 

Typically, to avoid the task of having to track/capture 
shocks in the stress/velocity field, approximations to 
MCP are invoked which give continuous solutions ei- 
ther by altering the boundary conditions or simplifying 
the PDE's. Smooth stress approximations are especially 
useful when attempting to solve for the velocity field — 
tracking flow shocks coming from both a discontinuous 
stress field and hyperbolicity in the velocity equations is 
an enormous job. To our knowledge, a full solution to 
MCP has never been obtained either numerically or an- 
alytically in cases where the underlying stress field has 
shocks. Instead, shock-free approximate solutions have 
mostly been pursued. 

Arguably, the two most successful and commonly used 
results of MCP are actually approximations, not full so- 
lutions. The Jenike Radial Solution |40|, |4^ for wedge 
hopper flow solves the MCP equations exactly, and with 



no discontinuities, but does not allow for a traction free 
top surface. It is a similarity solution of the form 

P = rf{e) (19) 
^ = 9{0) (20) 



which reduces the entire system to 3 ODE's with (r, 0) 
the position (r = distance from the hopper apex and 
is measured anticlockwise from the vertical). Though 
this solution enables the material to obey a wall yield 
criterion along the hopper walls, the stress parameters 
at the top surface have very little freedom. This is why 
most claim the Radial Solution to only hold near the 
oriflce, considerably away from the actual top surface. 

Another commonly used simpliflcation is called the 
Method of Differential Slices, although it only applies 
to stresses and not flow (our focus here). Originally pro- 
posed by Janssen in 1895 and signiflcantly enhanced since 
then, it is used to determine wall stresses in bins and 
containers. The method makes some very far-reaching 
assumptions about the internal stresses: p is presumed 
to only depend on height and the ip fleld is assumed to 
be identically 7r/2 or 0. These assumptions reduce the 
stress balance equations to one ODE and ultimately give 
the famous result that wall stresses increase up to a cer- 
tain depth and then saturate to a constant value. (This 
saturation behavior is not a byproduct of the approxima- 
tion; the discontinuous, full solution to the stress balance 
equations in a bin also gives similar stress saturation be- 
havior.) While this effect has been verifled extensively in 
experiment, the underlying assumptions clearly cannot 
hold since, for example, the walls exert an upward shear 
stress on the material which contradicts the assumption 
about tp. j^. 

In summary, the equations of MCP theory have very 
limited applicability to granular flows. There are very 
few, if any, solutions available (either numerical or an- 
alytical) for many important geometries such as planar 
or annular Couette cells, vertical chutes, inclined places, 
etc. In the case of silos, MCP has been used extensively 
to describe stresses, although the equations are difficult 
to solve and poorly behaved from a mathematical point 
of view, as noted above. There have also been some at- 
tempts to use MCP to describe granular drainage from 
silos, in conjunction with the coaxial fiow rule, but this 
approach has met with little success, as we now elabo- 
rate. 



B. Physical difficulties with coaxiality 

It is instructive to review the existing picture of silo 
drainage in MCP theory, to highlight what we will view 
as a major concern in the use of coaxiality for granular 
fiows. Suppose we have a fiat-bottomed quasi-2D silo 
with smooth side-walls. Under standard filling proce- 
dures, the walls provide only enough pressure to keep 
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(a) 



(b) 
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FIG. 8: Major principal stress chains in a quasi- 2D silo for 
the (a) active case, and (b) passive case. 



particles from sliding farther out. These wall conditions, 
known as the "active case" , give the following solution to 
the stress balance equations as found by marching down 
characteristics starting from the flat, pressure free, top 
surface: 



p{x,y) 



7t/2 

fgV 
1 + sin ( 



(22) 
(23) 



T : E 

-(Ji 

-0-3 

= \i\{<^3 - o-i) 

< 



l7l 
-|7| 



where A : B = T^AijBij and ±|7| are the principal 
strain-rate values. This type of behavior violates the 
second law of thermodynamics as it implies that mate- 
rial deformation does work on the system and likewise 
is non-dissipative. In more advanced plasticity theories, 
the thermodynamic inequality is upheld by requiring A in 
equation (P3|) to be non-negative, but in the basic limit- 
state framework we discuss, this contraint cannot be di- 
rectly enforced. 

We should quickly mention that in constructing the 
limit-state stress field for the discharging silo, we have 
used as a boundary condition that flow ensues when the 
pressure p above the hole drops differentially from the 
value it takes when the hole is closed. This claim allows 
us to preserve the continuous stress field described in 
equations and (|23|) for slow, quasi-static flow. 



where fg = pg is the weight density of the material and 
y is positive downward. Since the tjj field is identically 
7r/2 everywhere, the slip-lines are thus perfectly straight 
lines angled at ±e from the vertical. 

Refer again to Figure [SI The material deforms based 
solely on principal plane alignment. For a slow, dense 
flow in the silo geometry, coaxiality is troubling. Since 
the major principal stress is everywhere vertical, coaxi- 
ality requires material to stretch horizontally making it 
geometrically impossible for the material to converge and 
exit through the silo orifice. This issue has been dodged 
historically by asserting that a drastic change in the wall 
stresses occurs once the orifice opens, such that the walls 
drive the flow, not gravity 1]. The silo supposedly enters 
a "passive state" where the walls are pushing so hard 
on the material that the material is literally squeezed 
through the orifice by the walls. Even with this heavy- 
handed presumption, the solution predicted by equation 
(|T^ is absurd; it predicts the only non-stagnant regions 
in the silo are two narrow, straight channels which con- 
verge on the silo opening and are angled at ±45^ from 
the vertical. 

Coaxiality also suffers thermodynamically. The equa- 
tion itself only ensures there is no shear strain-rate in the 
principal stress reference frame and actually does not di- 
rectly enforce that of the two principal strain-rate axes, 
the axis of maximal compression (i.e. the major principal 
strain-rate direction) must align with the major princi- 
pal stress direction, as was the physical intuition shown 
in Figure [5J Coaxiality just as easily admits solutions 
for which the minor principal stresses align with the ma- 
jor principal strain-rate. When this happens, the plastic 
power dissipated per unit volume can be written 



C. Incipient yield everywhere 

The constraint of assuming a limit-state stress field of 
incipient yield everywhere is also questionable. Granular 
flows can contain regions below the yield criterion within 
which the plastic strain-rate must drop to zero. For ex- 
ample in drainage from a wide silo with a small orifice, 
the lower regions far from the orifice are completely stag- 
nant |43, |49] , and thus can hardly be considered to be 
at incipient yield. In fact, discrete-element simulations 
show that grains in this region essentially do not move 
from their initial, static packing ^2^. Simulations also 
reveal that high above the orifice, where the shear rate 
is reduced, the packing again becomes nearly rigid [Hoj . 
suggesting that the yield criterion is not met there either. 
As the silo example illustrates, a more general description 
of stresses coming from an elasto-plasticity theory may 
be necessary to properly account for sub-yield material 

Elasto-plasticity also alleviates another major concern 
with MCP which is that it is only well-defined in two 
dimensions. Three-dimensional stress tensors have six 
free variables, too many to be uniquely described by 
just force balance and incipient failure (altogether only 
4 equations). We mention in passing that extensions of 
MCP to axisymmetric three-dimensional situations have 
been developed. For example, the Har Von Karman hy- 
pothesis, which assumes that the intermediate principal 
stress (72, where ai < (J2 < era, is identical to either cti 
or (73, is frequently utilized in solving for conical hop- 
per flow. However, elasto-plasticity can handle three di- 
mensions without this hypothesis, while also allowing for 



9 




FIG. 9: A meso-scale object containing a small number of 
randomly packed, discrete grains, which controls the dynam- 
ics of dense granular flow, analogous to the "representative 
volume element" in classical continuum mechanics. 



stress states below the yield criterion in different regions. 



D. Neglect of discreteness and randomness 

Beyond its practical limitations and mathematical dif- 
ficulties, the most basic shortcomings of MCP are in its 
assumptions. Above all, a granular material is not con- 
tinuous. The microscopic grains composing it are usually 
visible to the naked eye, and significant variations in the 
velocity often occur across a distance of only several par- 
ticle lengths, e.g. in shear bands and boundary layers. 
Of course, the general theory of deterministic continuum 
mechanics is only expected to apply accurately when the 
system can be broken into "Representative Volume El- 
ements" (RVE's) of size L fulfilling d <^ L <^ Lmacro 
for d the microscale (particle size) and Lmacro the size of 
the system [H^ , which is clearly violated in many granu- 
lar flows. Therefore, the discrete, random nature of the 
particle packing must play an important role in the defor- 
mation process. To incorporate this notion coherently, it 
may be useful to seek out a dominant meso-scale object 
as a substitute for the RVE, upon which mechanical flow 
ideas apply, but in a non-deterministic, stochastic fashion 
(see Figure This concept is somewhat comparable to 
the "Stochastic Volume Element" in the theory of plas- 
ticity of heterogeneous materials 01 • that setting, it 
is known that (what physicists would call "quenched") 
randomness in material properties leads to blurring of 
the slip-lines, but, to our knowledge, this effect has not 
been considered in MCP for granular materials. 

More importantly, however, since the meso-scale 
should only be a few grains in width, there must also 
be randomness in the dynamics of yielding to an applied 
stress or body force, since the theoretical concept of a 
continuous slip-line is incompatible with the reality of 
a discrete, random packing. A stochastic response also 
seems fundamentally more consistent with the assump- 
tion of inicipient yield: If the material is just barely in 
equilibrium, it must be very sensitive to small, random 
fluctuations, causing localized yielding. 

We conclude that the shortcomings of MCP may be 



related to the deterministic Coaxial Flow Rule, so we now 
proceed to replace it with a more physically appropriate 
Stochastic Flow Rule. The Mohr-Coulomb description 
of stresses is more clearly grounded in principles of solid 
mechanics and is widely used in silo engineering, so we 
will assume that it still holds, on average, in the presence 
of slow flows, as a first approximation. 



IV. THE STOCHASTIC FLOW RULE 
A. Diffusing "spots" of plastic deformation 

It has been noted in a variety of experiments that dense 
granular flows can have velocity profiles which ressem- 
ble solutions to a diffusion equation. By far the best 
example is drainage from a wide silo, which has a well- 
known Gaussian profile near the orifice, spreading ver- 
tically as the square root of the height (with parabolic 
streamlines) in a range of experiments 26,^,^, 53, 
Recently, experiments in the split-bottom Couette geom- 
etry have demonstrated precise error-function profiles of 
the velocity spreading upward from the shearing circle, 
albeit with more complicated scaling 55]. Shear bands, 
when they exist, tend to be exponentially localized near 
moving rough walls, but we note that these too can be 
viewed as solutions to a steady drift-diffusion equation 
with drift directed toward the wall. 

It seems, therefore, that a successful flow rule for dense 
granular materials could be based on a stochastic process 
of deformation, consistent with our general arguments 
above based on discreteness and randomness. This begs 
the question: What is the diffusing carrier of plastic de- 
formation? In crystals, plasticity is carried by disloca- 
tions, but it is not clear that any such defects might exist 
in an amorphous material. The Gaussian velocity profile 
of granular draina ge w as first explained independently 
by Litwiniszyn [23. l56j and Mullins 25] as being due to 
the diffusion of voids upward from the orifice, exchanging 
with particles to cause downward motion. However, this 
model cannot be taken literally, since granular flows have 
nearly uniform density with essentially no voids and with 
far less cage-breaking than the model would predict. 

Instead, the starting point for our theory lies in the 
work of Bazant [^^l, who proposed a general model for 
the flow of amorphous materials (dense random packings) 
based on diffusing "spots" of cooperative relaxation, as 
illustrated in Figure^] The basic idea is that each ran- 
dom spot displacement causes a small block-like displace- 
ment of particles in the opposite direction. This flow 
mechanism takes into account the tendency of each par- 
ticle to move together with its cage of first neighbors, so 
the size of a spot is typically three to five particle diame- 
ters, Ls ^ 3 — 5d. This expectation has been confirmed in 
dense silo drainage as the length scale for spatial veloc- 
ity correlations in the experiments of Choi et al. [4^ . IH^ 
(data shown in Figure [TT|) as well as the discrete-element 
simulations of Rycroft et al. |2^]. In the continuum me- 
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FIG. 10: Spots as carriers of plastic deformation in amor- 
phous materials, (a) Cartoon of basic spot motion. A spot 
of local fluidization, carrying some free volume, moves to the 
upper right causes a cooperative displacement of particles, on 
average to the lower left, opposing the spot displacement, (b) 
In silo drainage, spots are injected at the orifice and perform 
random walks biased upward by gravity, causing downward 
motion of particles, (c) In other situations, we conjecture 
that spots are created during initial shear dilation and per- 
form random walks biased by local stress imbalances and body 
forces during steady flow. 




FIG. 11: Spatial velocity correlations in silo drainage exper- 
iments as in Refs. pfll . l57j with glass beads (d — 3mm) ob- 
tained high-speed digital-video particle tracking. The cor- 
relation coefficient of instantaneous displacements of a pair 
different particles is plotted as a function of their separation, 
averaged over all pairs and all times in the video. (Courtesy 
of the authors of Ref. 57].) 



spots), the spot simulations were able to accurately re- 
produce the statistical dynamics of several hundred thou- 
sand frictional, visco-elastic spheres in discrete-element 
simulations of drainage from a wide silo. This suggests 
that a general understanding of dense granular flows may 
come from a mechanical theory of spot dynamics. 



B. General form of the flow rule 



chanics terminology, this suggest that the spot may be 
an appropriate meso-scale replacement for the RVE. 

A major motivation for our work comes from the re- 
cent demonstration by Rycroft et al. that the Spot Model 
can be used as a basis for realistic multiscale simulations 
of dense granular drainage in a wide silo, assuming that 
spots perform upward random walks, biased uniformly 
by gravity 28]. Using only five fitting parameters (the 
size, volume, diffusivity, drift speed, and creation rate of 



In this work, we propose such a mechanical theory, 
based on the assumption that MCP provides a reasonable 
description of the mean stresses in slow dense granular 
flows. The key idea is to replace the Coaxial Flow Rule 
with a "Stochastic Flow Rule" based on mechanically 
biased spot diffusion. In the continuum approximation, 
the general form of the flow rule thus consists of two 
steps |2^: (i) a Fokker-Planck (drift-diffusion) equation 
is solved for the probability density (or concentration) of 
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spots, ps{r,t): 



dt ^ dx- 



_d d_ 



(of 



(24) 



where {D^} is the drift vector and {D^^} the diffusivity 
tensor of spots, determined by the stress field (below); 
and (ii) the mean drift velocity of particles u = {u"^} is 
constructed to oppose the local flux of spots: 



= - J dr'w{r,r')[D'^{r',t)ps{r',t) 



d 



[Df{r',t)ps{r',t)) 



(25) 



where w{t^ r') is the (dimensionless) spot "influence func- 
tion" specifying how much a particle at r moves in re- 
sponse to a spot displacement at r^ Without making 
the continuum approximation, the same physical picture 
can also be the basis for a multiscale model, alternating 
between macroscopic continuum stress computation and 
discrete spot-driven random-packing dynamics j28j . 

The mean flow proflle is derived from a nonlo- 
cal stochastic partial differential equation for spot-driven 
particle dynamics, in the approximation that spots do 
not interact with each other ^22}- Here, we have as- 
sumed the centered Stratonovich deflnition of stochastic 
differentials 58], which means physically that the spot 
influence is centered on its displacement. In contrast, 
Bazant used the one-sided reverse-Ito deflnition, where 
the spot influence is centered on the end of its displace- 
ment, which leads to an extra factor of two in the last 
term This choice is mathematically unrestricted 

(the "stochastic dilemna" |5^), but we flnd the centered 
deflnition to be a somewhat more reasonable physical hy- 
pothesis. Rycroft et al. have also found that the centered 
deflnition produces more realistic flowing random pack- 
ings in multiscale spot simulations of granular drainage, 
when compared to full discrete-element simulations 28]. 
If we use the simple approximation w ^ d{\T — t'\) in the 
integral for particle velocity, the Stratonovich interpreta- 
tion has the beneflt of automatically upholding volume 
conservation. 

Without even specifying how local stresses determine 
spot dynamics, the general form of the flow rule 
predicts continuous velocity flelds, even when the mean 
stresses are discontinuous. For example, shocks in the 
MCP stress fleld may lead to discontinuities in the spot 
drift vector, Di, and thus the spot flux. However, due 
to the nonlocal nature of the spot model, the particle 
flux is a convolution of the spot flux with the spot influ- 
ence function, thus preserving a continuous velocity fleld, 
which varies at scales larger than the spot size, Lg. This 
is a direct consequence of the geometry of dense random 
packings: The strong tendency for particles to move with 
their nearest neighbors smears velocity changes over at 
least one correlation length. 

In the simplest approximation, the spot influence is 
translationally invariant, w = w;(r — r'), so that spots 



everywhere in the system have the same size and shape. 
The spot influence decays off for r > L^, as a Gaussian 
among other possibilities. 

In we allow for the likelihood that the spot in- 
fluence may not be translationally and rotationally in- 
variant, e.g. since the local stress state always breaks 
symmetry. This is actually clear in the experimental 
measurements of Figure ^2 where velocity correlations 
are more short-ranged, without roughly half the decay 
length, in the vertical direction (parallel to gravity) com- 
pared to the horizontal (transverse) direction. This sug- 
gests that spots are generally non-spherical and may be 
more elongated in the directions transverse to their drift 
(or the body force). If anisotropy in the spot influence 
were taken into account, it would also be natural to allow 
for an anisotropic spot diffusivity tensor, which depends 
on the local stress state. However, such effects seem to 
be small in the granular flows we consider below, which 
are well described by a much simpler model. 

Another likely possibility is that spots come in a range 
of shapes and sizes. There could be a statistical distribu- 
tion of regions of local fluidization or plastic yield related 
to the local packing and stress state. It is thus more rea- 
sonable to think of the spot influence function w{r^r') as 
averaging over this distribution, just as a spatial velocity 
correlation measurement averages over a large number 
of collective relaxations. One advantage of taking the 
continuum limit of a Fokker-Planck equation in ap- 
plying the SFR is that such details are buried in the coef- 
flcients, which could in in principle be derived systemat- 
ically from any microscopic statistical model, or simply 
viewed as a starting point for further physical hypotheses 
(as we do below). 

Finally, we mention that there are also surely some 
nontrivial interactions between spots, which would make 
the SFR nonlinear and could lead to interesting new phe- 
nomena, such as spontaneous pattern formation. For ex- 
ample, spots may have a medium range attraction, since 
it is more diflicult to propagate particle rearrangements 
and plastic deformation into less dense, less mobile re- 
gions; there could also be a short range repulsion if the 
spot density gets too high, since grains will collapse into 
overly dilated regions. Such effects may be responsible 
for intermittent density waves in draining hoppers with 
narrow oriflces ^ , and perhaps even shear banding 
in other amorphous materials, such as metallic glasses, 
with a different plastic yield criterion. However, we will 
see that the hypothesis of non-interacting spots already 
works rather well in cases of steady, dense granular flows. 



C. A simple model for steady flows 

Due to efficient dissipation by friction, granular ma- 
terials subjected to steady forcing typically relax very 
quickly to a steady flowing state. For example, when 
a silo's oriflce is opened, a wave of reduced density 
(spots) progates upward, leaving in its wake a nearly 
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steady particle flow, which we associate with a steady 
flow of spots. This initial density wave can be seen very 
clearly in discrete-element simulations of various hopper- 
silo geometries 50j. For a narrow orifice, we have noted 
that intermittent flows with density waves can be ob- 
served bu t ty pical drainage flows are rather 
steady in time EJ • Similarly, when a rough inner 
cylinder is set into uniform rotation in a Couette cell, 
shear dilation propagates outward, raising the level of 
the packing, until a steady flow profile is reached. We 
interpret this initial dilation as signaling the creation of 
spots on the cylinder, which quickly reach a steady dis- 
tribution in the bulk. 

Hereafter, we focus on describing steady flow profiles, 
with equilibrium spot densities. For simplicity, let us as- 
sume isotropic spot diffusion, D'^^ = D2Saf3, since fluctu- 
ations are dominated by the (largely isotropic) geometry 
of dense random packings. Using the spot size Ls as the 
natural length scale, we can express the spot drift speed, 
|Di| = Lg/Ati, and diffusivity, D2 = L^/2At2, in terms 
of the times, Ati and At2, for drift and diffusion to reach 
this length. 

The flow profile of a draining silo, normalized by the 
outflow speed, is approximately constant over a wide 
range of flow speeds, as has recently been verified to 
great precision in the experiments of Choi et al. 57]. 
Not only is the mean velocity profile independent of flow 
rate (over an order of magnitude in mean velocity), but 
fluctuations about the mean, such as vertical and hor- 
izontal diffusion and measures of "cage breaking", also 
depend only on the distance dropped, and not explicitly 
on time (or some measure of "granular temperature"). In 
statistical terms, changing the flow rate is like watching 
the same movie at a different speed, so that the random 
packing goes through a similar sequence of geometrical 
configurations regardless of the velocity. Similar features 
have also been observed in shearing experiments in Cou- 
ette cells and numerical simulations of planar shear 

The experimental and simulational evidence, therefore, 
prompts the crucial approximation that Ati/ At2 = con- 
stant so as to uphold statistical invariance of the par- 
ticle trajectories under changing the overall flow rate. 
This can be justified if spots perform random walks 
with displacements selected from a fixed distribution, 
set by the geometry of the random packing ^2^. Here, 
we make the stronger assumption that the characteris- 
tic length of these random walks is the spot size L5, so 
that Ati = At2 = At. Our physical picture is that a 
spot represents a "cell" of localized fluidization (or plastic 
yield) of typical size L5, which triggers further fluidiza- 
tion ahead of it and randomly propagates to a neighbor- 
ing cell of similar extent. This picture is also consistent 
with the interpretation of Ls as a velocity correlation 
length above. 

With these hypotheses, the Fokker-Planck equation 



takes the simple time-independent form, 

where ds(r) = Di/|Di| is the spot drift direction, deter- 
mined by the mechanics of plastic yielding (below). The 
flow field is then 

u = 1 dr'w{r,r') (^d,(r')/C,(r') - ^Vp.(r')) 

(27) 

Equations (|26|) and ((?7)) define a simplified Stochastic 
Flow Rule, with only one parameter, L5, which need not 
be fitted to any flow profile. Instead, it can be measured 
independently as the velocity correlation length, which 
may be viewed as a dynamical material property. 



D. A mechanical theory of spot drift 

The main contribution of this paper is a simple the- 
ory connecting the spot drift direction to the stresses 
in MCP. The basic idea is to view the displacement of 
a spot as being due to a local event of material failure 
or fluidization. To make a quantitative prediction, we 
first define a cell of the material as the roughly diamond 
shaped region encompassed by two intersecting pairs of 
slip-lines, separated by Ls. When a spot passes through 
this cell, it fluidizes the material and thus locally changes 
the friction coefficient from the static value /i to the ki- 
netic value /i/e. This upsets the force balance on the cell 
and may cause a perturbative net force to occur. 

The force diagram for a material cell occupied by a 
spot can be broken into the sum of two diagrams (Figure 
[T^ . one which is the static diagram multiplied by /i/c/M 
and one which contains only normal contact force contri- 
butions and a body force term. MCP requires the static 
diagram to be balanced, thus the latter is the only cause 
for a net force. A well-known corollary of the divergence 
theorem enables us to express the surface integral of nor- 
mal contact stresses in terms of a gradient of p giving 



Fnet = (1 - ^) {^body - COS^ </>Vp) (28) 

as an effective force which pulls on a cell as it is fluidized 
by a passing spot, causing the spot to preferentially drift 
in the opposite direction. 

A spot cannot move in an arbitrary direction, however, 
since the material is at incipient yield only along the two 
slip lines. Therefore, the net force is constrained to pull 
the material cell along one of the slip-lines. The spot 
drift direction is then obtained by projecting (minus) the 
force, —Fneti onto the slip-lines and averaging these two 
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FIG. 12: (a) Material cell in static equilibrium. (b) A spot 
enters the cell fluidizing the material. In the force diagram, 
this means /i decreases to /Jk- (c) The force diagram for the 
fluidized material cell is best analyzed by breaking it into the 
sum of two diagrams. 



projection vectors with equal weight: 



(29) 
(30) 



where hg = (cos sin6>). With a formula for now 
determined, the SFR as stated in equations and (|77)) 
is now fully defined and ready for use. 

This continuum mechanical theory of spot drift also 
helps us understand the sources of spot diffusion. As 
noted above and sketched in Figure El a material cell is 
a small fragment of a random packing, which is unlikely 
to be able to accomodate shear strain precisely along the 



slip-lines of the mean continuum stress field. Instead, 
the instantaneous slip-lines are effectively blurred by the 
discrete random packing. Still, we preserve the picture 
of spots moving along slip-lines in constructing dg, but 
represent the additonal bluriness in the slip-line field by 
enforcing isotropic spot diffusivity. 



E. 



Frame indifference 



Finally, we must check that our flow rule satisfies frame 
indifference; solving for flow in different rigidly moving 
reference frames cannot give different answers in a fixed 
reference frame. Since the SFR is a 2D steady-state flow 
rule, the only flows we need to check for indifference are 
those with rotational/translational symmetry. In these 
cases, the particle velocity is a function of only one spa- 
tial variable and equation (|26|) for pg becomes a second- 
order ODE. In solving the boundary value problem, we 
must ensure that grains along the walls have a velocity 
vector tangent to the walls. This constrains one of the 
two degrees of freedom in the set of possible solutions 
for ps. Since is homogeneous, the other degree of 
freedom must come out as a multiplicative undetermined 
constant. Thus the velocity profile is unique up to a mul- 
tiplicative constant. 

With only one constant, we cannot match boundary 
conditions for particle speed along more than one wall in 
general. So to solve for a flow between two walls, we must 
add rigid-body motions to the reference frame of the ob- 
server until we have the unique frame for which a solution 
exists matching both boundary conditions. This is an un- 
expected and welcome bonus of the SFR. Most flow rules 
in continuum mechanics enforce material frame indiffer- 
ence directly, i.e. the flow rule itself is derived to be au- 
tomatically satisfied by any rigid-body motion, ensuring 
the same results independent of reference frame. Coax- 
iality achieves this by relating stress information only 
to strain-rate variables for instance. The SFR, however, 
upholds frame indifference indirectly in that the solution 
does not exist unless the problem is solved in exactly one 
"correct" frame of reference. 

We have thus integrated the spot concept with the the- 
ory of plastic stresses treating spots as the "carriers of 
plasticity" . We note that up to our granular-specific de- 
termination of the drift direction, the SFR principle can 
be applied to any amorphous isotropic material with a 
small characteristic length scale (dominant randomness) 
and a yield criterion. 



V. APPLICATIONS TO GRANULAR FLOW 

The Stochastic Flow Rule is quite general and in prin- 
ciple can be applied to any limit-state plasticity model 
of stresses, with different choices of the yield function to 
describe different materials. In this section, we apply the 
simplest SFR ^^-^7^ to granular materials with MCP 
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stresses and compare its predictions to a wide range of 
existing experimental data for steady dense flows. In 
calculating stresses, we assume a typical friction angle 
of ^ = 30^. It is known that for spherical grains, the 
friction angle usually lies in a somewhat narrow range 
of about 20^ — 30^ and can be as lame as 50^ for some 
anisotropic, highly angular materials 1]. In the examples 
we consider, however, varying the (j) value in this range 
has very little macroscopic effect in our model. 

The spot size Lg has a much larger effect, so we fo- 
cus on its role in a variety of dense flows. We emphasize 
that we do not fit Lg to any flow profile below. Instead, 
we simply use the range Ls = 3 — 5(i for dense flowing 
sphere packings inferred independently from spatial ve- 
locity correlations in silo drainage experiments 57] and 
simulations 28] (see Fig. [TTj) . This is consistent with 
our view of the correlation length, Lg, as a fundamental 
geometrical property of a flowing granular material. 

Without any fitting parameters, we will apply the sim- 
ple SFR to several prototypical flows. Each has different 
forcing and symmetries and, to our knowledge, they can- 
not be simultaneously described by any existing model. 
Our first example is granular drainage to a small orifice 
in a wide fiat-bottomed silo, driven entirely by gravity. 
Our second example is shear flow in an annular Cou- 
ette cell driven by a moving rough inner cylinder, where 
gravity plays no role. Our third example is the drag- 
ging of a loaded plate over a semi-infinite material at 
rest, which combines gravitational forces and boundary 
forcing. Lastly we apply the SFR to a canonical free- 
surface flow, the continuous avalanching of a granular 
heap. The transition to a rapidly flowing surface shear 
layer on a heap will also lead us into a discussion of how 
rate-dependent effects, such as Bagnold rheology, may 
naturally extend into our model. 

Throughout our treatment of the various examples, the 
first step will be to solve (|26|) to obtain the "unconvolved" 
velocity field 



l2 



(31) 



which corresponds to the SFR velocity ((Tf)) for a point- 
like influence function w = 5{\r — r'\). For the most part, 
u* is the "skeleton" for the full solution u because con- 
volving u* with a general spot influence merely blurs out 
the sharper features of u*. In some situations with wide 
shear zones, such as silo flow, the convolution has only a 
minor effect, but in others with narrow shear bands, at 
the scale of the spot size, the convolution step is essential 
for self-consistency and accuracy. 



A. Silos 

The flow proflle in a flat-bottom silo geometry is well- 
known for its striking similarity to the fundamental so- 
lution of the diffusion equation. As noted above, early 
models of silo flow explained this based on the upward 



diffusion of voids from the oriflce [2J, [23, • Without 
reference to a speciflc microscopic mechanism, Nedder- 
man and Tiizun later derived the same equations based 
on a continuum constitutive law 11^. They asserted 
that the horizontal velocity component u is proportional 
to the horizontal gradient of the downward velocity com- 
ponent V, 



u = b 



dv 
dx 



(32) 



since particles should drift from regions of slow, dense 
flow toward regions of faster, less dense (more dilated) 
flow. Assuming small density fluctuations, mass conser- 
vation applied to the 2D velocity fleld, u = (i^, —v) then 
yields the diffusion equation for the downward velocity. 



dv ^ d'^v 
dz dx'^ 



(33) 



where the vertical direction z acts like "time". The dif- 
fusivity b is thus really a "diffusion length" , to be deter- 
mined empirically. An advantage of the continuum for- 
mulation is that it avoids the paradox (resolved by the 
Spot Model 27]) that the classical picture of void ran- 
dom walks requires 6 <C (i, while experiments invariably 
show b > d. 

Solving the Kinematic Model in the wide flat-bottomed 
silo geometry with a point oriflce gives the familiar Green 
function for the diffusion equation. 



v{x^ z) 



^-x^ /Ahz 

V4:7rbz 



(34) 



This gives an excellent match to experimental data close 
to the oriflce (e.g. see Figure IT^ . although the flt grad- 
ually gets worse with increasing height, as the flow be- 
comes somewhat more plug-like. Nevertheless, Gaussian 
flts of experimental data have provided similar estimates 
of 6 = 1.3d 57J, 1.3 - 2.3d 49], 2.3d 26], 3.5d 48], and 
2d — 4:d for a variety of granular materials composed 
of monodisperse spheres. 

We now apply our theory to this flow geometry and 
see how it connects to the Kinematic Model. Applying 
equation (pUj) using the stress fleld described by equations 
and gives uniform upward spot drift; Fnet comes 
out as pointing uniformly downward and the slip-lines are 
symmetric about the vertical (see Figure IT^ . The SFR 
(|26j) then reduces to 



dps 
dz 



2 



dh 



^ d'ps 

dx'^ dz'^ 



(35) 



although we emphasize that this form applies only when 
the walls are smooth or equivalently when the silo width 
is large. The last term, which represents vertical diffu- 
sion of spots (relative to the mean upward drift), makes 
this equation for the spot density differ somewhat from 
the simple diffusion equation for the downward velocity 
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FIG. 13: The mean velocity profile in a wide quasi-2d silo of 3mm glass beads from Ref. |4^. Horizontal slices of the downward 
velocity component near the orifice, indicated in the complete flow profile on the left, are shown on the right, and compared 
to the predictions of the Kinematic Model with two choices of the parameter b. The Stochastic Flow Rule for MCP for a wide 
silo (without side walls) gives a similar velocity field with 6 ~ 1.5 — 2.5d. 




Outflow 



FIG. 14: The flat-bottomed silo geometry. The intersect- 
ing black lines represent the slip-line field as determined from 
solving the stress balance equations of MCP, and the red vec- 
tor field is the spot drift direction, as determined from the 
SFR. In this highly symmetric geometry, the spot drift pre- 
cisely opposes the gravitational body force, d = —g. 



of the Kinematic Model. Consistent with our model, ver- 
tical diffusion, with a similar (but not identical) diffusion 
length as horizontal diffusion, has been observed in recent 
silo-drainage experiments 49, 57]. 

The general solution of can be expressed as a 
Fourier integral. 



1 

2^ 



Akx 



where A{k) is the Fourier transform of the spot density 
at the bottom {z = 0). The narrowest possible orifice al- 
lowing for flow is the case of a point source of spots, 
^5(^,0) (X S{x), ^{k) oc 1 (which is also the Green 
function). Convolution with a spot influence function 
of width Ls produces a downward velocity profile on the 
orifice of width Lg. Unlike the Kinematic Model (or any 
other continuum model which does not account for the 
finite grain size), our theory thus predicts that flow can- 
not occur unless the orifice is at least as wide as one spot, 
Ls = 3 — 5d. 

The details of flow very close to the orifice, z = 0{Ls), 
are controlled by the choice of boundary condition, re- 
flecting the dynamics of dilation, contact-breaking, and 
acceleration at the orifice, which are not described by 
our bulk dense-flow model. Rather than speculate on the 
form of this boundary condition, we focus on the bulk re- 
gion slightly farther from the orifice. For z ^ Lg (and 
Lsk <C 1), the vertical diffusion term becomes unimpor- 
tant, and the Green function tends to a Gaussian 



v{x^ z) 



where the variance is 



g-xV2<7^(z) 



0{Ll). 



(37) 



(38) 



dk 



(36) 



(The second term is an offset from convolution with the 
spot influence function, which also depends on the choice 
of boundary conditions.) 

There has been no prior theoretical prediction of the 
kinematic parameter 6, which we interpret as the spot 
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diffusion length j^^. Comparing ((51)) and ((57)) . we ob- 
tain b = Ls/2 = 1.5 — 2.5d without any fitting, beyond 
the independent determination of Lg from velocity cor- 
relations. This prediction is in excellent agreement with 
the experimental measurements listed above. However, 
the model does not predict the apparent increase of b 
with height, as the fiow becomes more plug like. This 
may be due to the breakdown of the assumption of in- 
cipient yield higher in the silo, where the shear is greatly 
reduced, and it may require modeling stresses more gen- 
erally with elasto-plasticity. 

In any case, we are not aware of any other model of 
silo fiow with a plausible basis in mechanics. It is note- 
worthy that we assume active silo stresses (driven by 
gravity), as typically assumed in a quasi-static silo. As 
a result, we do not require a sudden switch to passive 
stresses (driven by the side walls) upon fiow initiation, 
as in existing plasticity theories based on the Coaxial 
Flow Rule jy. Our use of the standard MCP model for 
stresses in quasi-static silos also suggests that the SFR 
may predict reasonable dependences on the geometry of 
the silo or hopper, wall roughness, and other mechanical 
parameters. In contrast, the Kinematic Model fails to 
incorporate any mechanics, and, not surprisingly, fails to 
describe fiows in different silo/hopper geometries in ex- 
periments 49] . Testing our model in the same way would 
be an interesting direction for future work, since it has 
essentially no adjustable parameters. 

B. Couette cells 

The key benefit of our model is versatility; we will now 
take exactly the same model, which is able to describe 
wide silo fiows driven by gravity, and apply it to Taylor- 
Couette shear fiows in annular cells driven by a moving 
boundary. The granular material is confined between 
vertical rough-walled concentric cylinders and set into 
motion by rotating the inner cylinder. The fiow field has 
been characterized extensively in experiments and simu- 
lations, and several theories have been proposed for this 
particular geometry |23, IH^, • For example, a good 
fit of experimental data for Couette cells can be obtained 
by postulating a density and temperature dependent vis- 
cosity in a fiuid mechanical theory 20] , but it is not clear 
that the same model can describe any other geometries, 
such as silos, hoppers, or other shear fiows. 

To solve for the MCP stresses in the annular Couette 
geometry, we first convert the stress balance equations 
to polar coordinates (r, 6) and require that p and i/j obey 
radial symmetry. This gives the following pair of ODEs: 

or ^ sin2V>* 
dr r(cos 2ip^ + sin cj)) 

^ = (40) 

dr r(cos 2?/^* + sin (j)) 

where rj = \ogp and ?/^* = iIj -\- ^ — 0. Although ?/^* has 
an implicit analytical solution, r] does not, so we solve 



(a) 




FIG. 15: (a) A plan view of the annular Couette cell geometry, 
where the granular material is confined between concentric 
vertical cylinders. The rough wall is rotated ant i- clockwise 
while the outer wall is held fixed. The crossing black lines 
within the material are the slip-lines as found from MCP, 
and the red vector field is the spot drift as determined by 
the SFR. (b) Normalized SFR velocity from as a function of 
distance from the inner wall with inner cylinder radius 15d, 
25d, 50d, and lOOd (from bottom to top curves, respectively). 
The friction angle is = 30°, and the spot size is Ls = Sd. 



these equations numerically using fully rough inner wall 
boundary conditions '0*(r^(;) = f — ^ and any arbitrary 
value for rjivy^). The resulting slip-lines are shown in 
Figure irKT a). 

In the Couette geometry, the average normal stress, p, 
decreases with radial distance, which implies that the fiu- 
idization force on material, Fneti is everywhere directed 
outward. We then apply equation ((30)) to calculate the 
drift direction ds{r) by projecting the vector Fnet onto 
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slip-lines, and implement the SFR, exploiting symmetry 
which allows only a nonzero velocity in the direction. 
This implies 

u*.f = = -L,(d,-f)p, + :^^ (41) 

which yields a solution for pg up to a scalar factor. We 
then use ps to calculate the 9 component of the (uncon- 
volved) velocity once again using the SFR equation, 

u* J = -L,(d, -^jp,. (42) 

It turns out, as we may have expected, that u* has a 
shear band at the inner wall with nearly exact exponen- 
tial decay. The length scale of this decay is the spot size, 
Lg, since this is the velocity correlation length, beyond 
which the inadmissible shear at the inner cylinder can be 
effectively dissipated by the material. 

The thinness of the shear band requires that, for con- 
sistency, we must take into account the finite spot size in 
reconstructing the velocity field through the convolution 
integral ((77|) . For simplicity we will use a uniform spot 
influence function, i.e. 

where H{x) is the Heaviside step function. To evaluate 
the integral ((23), we also must make a hypothesis about 
how spots operate when they overlap one of the walls. 
Random packing dynamics near walls is different than in 
the bulk and sensitive to surface roughness, and further 
detailed analysis of experiments and simulations will be 
required to elucidate the collective mechanism(s). Here, 
the precise shape of spots near the wall has little effect, 
except to flatten out the spike in velocity that occurs 
near the wall in the unconvolved velocity u*. As sim- 
ple first approximation, used hereafter in this paper, we 
will view the space beyond each boundary as containing 
a bath of non-diffusive spots at uniform concentration 
whose flux is such that the particle velocity invoked "in- 
side" the boundaries directly mimics the rigid motion of 
the walls. This effectively "folds" part of the spot in- 
fluence back into the granular material when it overlaps 
with the wall. The resulting velocity field is shown in 
Figure ITM b). where normalized velocity is shown versus 
distance from the inner cylinder wall for Ls = 3d for a 
wide range of inner cylinder radii. 

The predicted flow field - without any fitting - is in 
striking agreement with a large set of data from ex- 
perimental and discrete-element simulations for different 
cylinder radii and grain sizes 19, 33, 63,^]. As shown 
in Fi gure [TBI the experimental data compiled by GDR 
Midi falls almost entirely within the predicted SFR 
velocity profiles, by setting the spot size to the same typ- 
ical range of correlation lengths, = 3 — 5(i, measured 
independently in a quasi-2D silo (Figure [TT|) . Viewing 
the data on a semilog plot shows that the agreement ex- 
tends all the way into the tail of the velocity distribu- 
tion. We emphasize that the same simple theory, with 



(a) 




FIG. 16: Theory versus experiment for the normalized ve- 
locity in annular Couette cells on (a) linear and (b) semilog 
plots. The dashed curve is the predicted SFR velocity field 
with Ls = 3d, while the solid line is for Ls = 5d; both curves 
are for an inner cylinder radius of rw = SOd and = 30°. 
Experimental measurements (points) for a wide range of in- 
ner and outer cylinder radii are shown from the compilation 
of data shown in Figure 3c of 33] . (The experimental data is 
courtesy of GDR Midi and originates from the work of [gJ, 
m, 11^, and £7J.) 



the same range of Ls values, also accurately predicts silo 
flows above, as well as other situations below. Unify- 
ing all of this data in a single simple theory without any 
empirical fitting constitutes a stringent quantitative test. 

It is interesting to note the behavior close to the wall, 
especially in thin Couette cells. In experiments [s^, an- 
nular flow profiles are known to have a Gaussian correc- 
tion term when the thickness of the cell becomes non- 
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FIG. 17: Predicted variation in the width of the shear band 
with SFR over the standard range of granular friction angles 
for the annular Couette geometry. (Lg = 3d, Vw — SOd) 



negligible in terms of particle size. This slight flatten- 
ing near the wall is apparent in our model as well and 
is a byproduct of convolving with the spot influence. 
We thus interpret this feature as another sign of the 
strongly correlated motion of particles, primarily with 
the "cage" of nearest neighbors, as approximately de- 
scribed by the spot mechanism. In this calculation, we 
used a uniform spot influence, but have noticed relatively 
little sensitivity of the predicted flow proflle, for different 
strongly localized influence functions, such as a Gaus- 
sian, w (X e~^^ A detailed comparison of the model 
to experimental data may provide fundamental insights 
into the spot influence, and thus the collective dynam- 
ics of random packings, near a rough wall at the discrete 
particle level. 

The experimental results shown in Figure^Jcome from 
apparati with inner wall radii ranging from lAd — lOOd. 
The relatively small variations in the data sets over such 
a large range of inner radii clearly indicates that the in- 
ner wall radius is not a crucial length scale in the flow. 
The plotted theoretical prediction uses an inner radius of 
~ 80(i, but, as can be seen in Figure ITCT b). our results 
depend only minimally on the inner cylinder radius. In- 
deed, the meso-scale correlation length of = 3 — 5(i is 
the dominant length scale in our theory for this geometry. 

To substantiate an earlier claim, we now consider how 
the friction angle (j) affects the flow properties (holding 
Ls flxed) according to our model. We can see this most 
clearly by observing how the shear band half-width (i.e. 
the distance from the wall to the location where velocity 
is half-maximum) varies over the (j) range for usual gran- 
ular materials (~ 20^ — 50^). As shown in Figure [T7I the 
half- width changes by < {)Ad over the entire range and by 
< 0.1(i for the range of laboratory-style spherical grains. 
This very weak influence of internal friction agrees with 
simulations in the Couette geometry by Schollmann |6^ . 



C. Plate dragging 

We now examine perhaps the simplest situation where 
gravity affects the shear band caused by a moving rough 
wall. Consider slowly dragging a rough plate horizontally 
across the upper surface of a deep (semi-inflnite) granular 
material. The plate maintains full contact by pressing 
down on the surface with pressure po cos^ (j). The proflle 
of the shear band which forms below the plate depends 
on the relative loading pressure, = Pq/ fgi where fg 
is the weight (gravitational body force) density of the 
material. 

The plate-dragging flow fleld can be found using a pro- 
cedure analagous to the annular Couette cell, but en- 
forcing horizontal instead of radial symmetry. With y 
measuring distance below the plate, the stress balance 
equations give 

— sin 2ip 
^ 2(7(cos2?/; - sin^) 
cos 2?/^ 
cos 2?/^ — sin (j) 

where q{y) = p{y)/fg is the average normal stress scaled 
to the weight density. The fluidization force will push 
material downward and spots upward resulting in a flow 
proflle that decays close to exponentially near the moving 
wall. 

Experiments jTol . FH and simulations [f^ offer 
differing assessments of the details of the flow proflle away 
from the shear band, but the dominant exponential de- 
cay behavior is clearly observed in all. The displayed SFR 
prediction (Figure [TH|) uses loading parameters from Tsai 
and Gollub y9] in order to appropriately compare with 
their results. Although the general properties of the flow 
appear to be represented well by the model, we do notice 
that the predicted range of typical flows is too small to 
fully encompass the experimental data. There could be 
a number of reasons for this discrepancy, but it is worth 
pointing out that the quasi-2D plate-dragging geometry 
is rather difficult to realize in experiments. For exam- 
ple, this experiment was executed by rotating a loaded 
washer-type object on top of an annular channel, and 
it was observed that the sidewalls pushing in the third 
dimension actually did play some role. 

In Refs. 73] and 75], simulations of this environment 
indicate that the shear band width increases with increas- 
ing loading of the top plate. As can be seen in Figure 
[T^ our theory captures this general trend of increasing 
loading causing increasing shear band width. However, 
the swing in band size predicted by our theory is not 
large enough to match the range of band sizes in simu- 
lations 73] and 75] in which the shear band half-width 
can be as large as several tens of particle diameters for 
large enough qq and diverges as go ^ oo (i.e. zero grav- 
ity) . In cases such as these where the value of qo becomes 
very large, as we will discuss in more depth after the next 
section, we believe a new phenomenon begins to domi- 
nate our meso-scale argument and that this phenomenon 
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FIG. 18: (a) The plate- dragging geometry. The top wall drags along the top of a bed of granular material. The crossing black 
lines within the material are the slip-lines as found from MCP, and the vector field is the spot drift as determined by the SFR. 

(b) Theory against experiment for the plate dragging geometry. Theory: L=3 ( ), L=5 ( — ). Experiment (*) courtesy of 

the authors of 69] . 



may be attributed to a particular property of the slip-line 
field. 



D. Slow heap flows 

We now examine a prototypical free surface flow. Very 
close to the respose angle, a granular heap which is slowly 
but consistently re-fed grains undergoes a particular type 
of motion characterized by avalanching at the top sur- 
face and a slower "creeping" motion beneath. This type 
of flow has been studied in experiments ^] and simu- 
lations • Though heap flows with faster top shear 
layers have also been studied HE Iz3 we will focus for 
now on the slower regime, which more closely resembles 
a quasi-static flow where the SFR might apply. 

This kind of flow is stable, but indeed quite "delicate" 
in the sense that relatively small changes to the system 
parameters (i.e. flow rate, height of the flowing layer) 
can invoke large changes to the qualitative flow profile 
especially in the top layers [73 • We will describe and 
attempt to explain this effect more in the next section. 

The heap geometry is depicted in Figure [501 along with 



the corresponding spot drift field and slip-line field. Any 
gravity driven free surface flow problem for which the 
stresses and flow are approximately invariant in the direc- 
tion parallel to the surface will have limit-state stresses 
that obey the following relations: 



^ = 
P = 



—e 

fgV 
cos 6 



(44) 
(45) 



where y is the depth measured orthogonally from the 
free surface. Note that in limit-state theory, for self- 
consistency, the static angle of repose is identical to the 
internal friction angle (/>, which is a reasonable assertion 
but still debated in the community. (By "static repose 
angle" we refer to the angle of inclination below which a 
flowing system jams; in simulations of flow down a rough 
inclined plane, it has been shown that this angle does 
vary in a narrow range depending on the height of the 
flow [zJI?^.) Applying equation (|30j) to these equations 
gives a simple expression for the spot drift vector: 
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(46) 
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FIG. 19: Plot of theoretical shear band half- width vs relative 
loading pressure — p{y — 0) / fg of the top plate (where fg 
is the weight density). The calculation assumes Lg — 5d 
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FIG. 20: (a) The heap flow setup, (b) The dashed rectangle 
in (a) is enlarged; the slip-line field from MCP is plotted along 
with the drift field from the SFR. 



We may then apply the SFR, which simplifies upon re- 
quiring that the flow run parallel to the free surface (i.e. 
u* = {-u,0)). 

Since the drift field is uniform, we obtain an analytical 
solution for the unconvolved velocity: 

Thus our model predicts that the velocity decays expo- 
nentially off the free surface. In cases like these, where 
the boundary of the flow makes no contact with a rigid 
wall, it is less clear how the spots (and free volume) might 
behave near the flowing free surface. To avoid addressing 
this issue in detail, we neglect convolution with the spot 
influence function and simply assume u ^ u*. 

In their experiments on slow heap flow, Lemieux and 
Durian have shown that the velocity profile in the 
flowing top layers is indeed well approximated by an ex- 
ponential decay. Furthermore, they found the flow in this 
regime to be continuous and stable. The decay law they 
obtained is 

V40)«exp(-^) 
which is very close to our predicted solution for Ls = 3d: 



' ^ ' 4.58(i' 

Silbert et. al report finding a similar decay profile 
at low flow rates in simulations of flow down a rough 
inclined plane, although the avalanching at the surface 
was intermittent. In conclusion, we have demonstrated 
a fourth, qualitatively different situation where the same 
simple MCP/SFR theory predicts the flow profile, with- 
out adjusting any parameters. 

VI. TRANSITION FROM THE SFR TO 
BAGNOLD RHEOLOGY 

A. Breakdown of the SFR 

In the last two examples, plate dragging and slow heap 
flow, there are limits where the SFR fails to predict the 
experimental flow profiles. In this section, we will ex- 
plain why the breakdown of the SFR is to be expected 
in these cases and others, whenever slip-lines approach 
"admissibility" and coincide with shear planes. In this 
singular limit of the SFR, we postulate a transition to 
Bagnold rheology. The stochastic spot-based mechanism 
for plastic yielding is thus replaced by a different physical 
mechanism, the free sliding along shear planes. 

For example, consider the case of plate dragging above. 
At large relative loading, the flow field resembles that 
of a zero-gravity horizontal shear cell (between shearing 
flat plates), and it appears that the SFR breaks down: 
With body forces and Vp both going to 0, equation ((Sn|) 
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gives Fnet = implying that spots have no drift and 
consequently the only SFR solution is u = 0. 

Problems also occur with flow down a rough inclined 
plane: Slightly increasing the flow rate (and consequently 
the flow height) or inclination angle causes the velocity vs 
depth relationship to exit the exponential decay regime 
detailed above and undergo significant changes, passing 
first through a regime of linear dependence 80, SI] to 
a regime resembling a 3/2 power-law of depth [TOL 183. 
IssI, opposite in concavity to the exponential decay 
regime. 

Why does the velocity profile for inclined plane flow un- 
dergo many different qualitative regimes depending del- 
icately on system parameters, while others (e.g. silos, 
annular cells) appear to be only weakly affected and al- 
most always exhibit the same (normalized) velocity pro- 
file? Tall inclined plane flows and zero-gravity planar 
shear flows have been successfully described in multiple 
experiments and simulations 62, 79, 82, 83, 85] by the 
empirical scaling law of Bagnold 16]. In this section, 
we suggest a means to reconcile and perhaps eventually 
combine these theories into a coherent whole. 



been employed as test cases for Bagnold scaling. In the 
planar-shear environment, the shear and normal stresses 
acting on the shear planes are spatially constant through- 
out the flow. Equation ((I^ therefore implies that the 
strain-rate is uniform; as a result, the velocity varies lin- 
early from one wall to the other. This result is known 
as Uniform Shear Flow (USF) and is easily verified in 
simulations of Lees-Edwards boundary conditions. For 
example, the rheology ((I^ has been demonstrated in the 
simulations of |85^. 

Applying Bagnold scaling to the inclined plane geom- 
etry, slightly above static repose, gives a shear stress ex- 
cess which grows linearly with depth and thus a shearing 
rate that grows as the square root of depth. This implies 
a velocity profile of the form 

U OC /l3/2 _ y3/2 (50) 

for y the depth variable and h the height of the flowing 
material (with no-slip bottom boundary condition). In 
this way, Bagnold scaling successfully explains the 3/2 
power law dependence noted above. 



B. Bagnold rheology 

Let us briefly review Bagnold's classical theory of gran- 
ular shear flow. In its original form, "Bagnold scaling" 
expresses a particular rate-dependency for granular flow 
whenever the solid fraction is uniform throughout: 

r(X7^ (48) 

where 7 is the rate of simple shear. To account as well for 
static stresses arising from the internal friction, a related 
variant of this scaling law is commonly used [Sq : 

r — /icr (X 7^. (49) 

It is in some sense a law for how the yield criterion can 
be exceeded when non-negligible strain-rates can absorb 
the extra shear stress. This constitutive law alone is an 
incomplete flow theory since it provides no way of pre- 
dicting whether or not the solid fraction will be uniform 
during flow or how a non-uniform solid fraction affects 
the above rheology. Bagnold originally explained the 
quadratic relationship between stress and strain-rate in 
terms of binary collisions as the joint effect of both the 
particle collision rate and the momentum loss per colli- 
sion being directly proportional to the strain-rate [T^ . 
Despite this collision-based argument, however, Bagnold 
scaling has been observed to hold well into the dense 
regime, whenever the solid fraction is approximately con- 
stant throughout the system. This seemingly contradic- 
tory observation can be justified in the hard-sphere limit 
(without body forces) by a Newtonian invariance argu- 
ment IH^, although it calls into question the underlying 
physical mechanism. 

Zero-gravity planar shear flow and thick inclined plane 
flow both exhibit nearly uniform density and thus have 



C. Slip-line admissibility 

The seemingly disparate flow mechanisms of the SFR 
and Bagnold rheology can be reconciled very naturally 
by considering the geometry of the slip-lines. In plastic- 
ity theory, all flows can be classified based on "slip-line 
admissibility". For admissible slip-lines, boundary con- 
ditions are such that the flow can, and presumably does, 
take place by continuous shearing along only one family 
of slip-lines. In mathematical terms, the slip-lines are ad- 
missible for a given flow, whenever the double-shearing 
continuum flow-rule ([T^ allows multiple solutions to the 
boundary value problem. 

Slip-line admissibility is the exception, not the rule, 
since it is highly unlikely that the prescribed velocity 
boundary conditions are fulfilled by a continuous shear 
on either slip-line family. Be that as it may, it so happens 
that planar shear flow and inclined plane flow are both 
slip-line admissible. This special property is shared by 
no other flow geometry studied in this paper, or, to our 
knowledge, elsewhere in the granular materials commu- 
nity. (Contrast the slip lines in Figure HTl with those in 
Figures [H and [ini) 

There is also an interesting difference in the density 
distributions. For admissible flows, the volume fraction 
is nearly uniform, and Bagnold rheology has a reason- 
able physical justification. For the more common case 
of inadmissible flows, as in silos and Couette cells, the 
volume fraction is typically highly nonuniform. In such 
cases, the SFR seems to provide an excellent description 
of the flow, and Bagnold rheology clearly does not apply. 

These observations motivate us to think of admissi- 
bility as a criterion for two very different microscopic 
mechanisms for granular flow: In admissible flows, ma- 
terial motion is a viscous dragging of material "slabs" 
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FIG. 21: (a) Slip- line field for gravity- free planar shear flow, 
(b) Slip-line field for inclined plane flow. Note that in both 
cases the shear planes will be aligned completely with one 
slip- line family. 



along one slip-line family (Bagnold dominated), whereas 
in inadmissible flows there is no clear choice as to which 
slip-line family should control the motion and thus ma- 
terial randomly chooses between both slip-line families 
(SFR dominated). Perhaps admissibility in the slip-line 
field causes the solid fraction to remain roughly uniform 
because the material in no sense has to collide head-on 
into neighboring material for it to move. In flows where 
non-uniform dilation does occur, experiments have shown 
the motion is nearly independent of any local Bagnold 
rheology, encouraging our strong distinction between the 
dynamics of flow problems of differing admissibility sta- 
tus 65, 87]. 

These considerations all lead us to the fundamental 
conjecture: 

Slip-line admissibility is a geometrical and 
mechanical indicator as to the relative im- 
portance of rate-dependency (Bagnold rheol- 
ogy) over rate-independency (SFR) in a dense 
granular flow. 

This means that a flow which has an admissible limit- 
state stress field will be dominated by rate-dependent 
effects when the yield criterion is only slightly exceeded. 



terms; slip-lines are defined by the quasi- static stresses 
as lines along which t — jia = 0, whereas shear planes 
are defined entirely by the velocity profile. In an admis- 
sible system, the shear planes coincide with one slip-line 
family. In inadmissible systems, the shear planes almost 
everywhere do not coincide with slip-lines. 

With admissible slip-lines, excess shear stress tends to 
be uniformly distributed throughout the system, result- 
ing in global Bagnold rheology. For example, consider a 
zero-gravity planar shear cell. If we were to apply ad- 
ditional shear stress to the body in a manner aligned 
with the admissible slip-line family, e.g. by increasing 
the wall shear stress above yield by some amount Ar, 
that additional shear stress would distribute itself within 
the material precisely along the slip-lines. Every horizon- 
tal slip-line within the bulk would thus receive a boost 
in shear stress of size Ar. In limit-state theory the slip- 
lines have the highest possible t — jia value a quasi-static 
material element can take - zero. Adding Ar additonal 
shear stress to the slip-lines means that r will maximally 
exceed jicj precisely along the slip-lines, and, by admissi- 
bility, every shearing plane. As a result, there will be a 
Bagnold contribution everywhere. Similarly, if we took a 
limit-state inclined plane geometry and increased the tilt 
angle some amount, an analagous boost in shear stress 
along the admissible slip-line family would occur causing 
r to exceed jicj precisely along all the shear planes. 

In contrast, with inadmissible slip-lines, excess shear 
stress tends to remain localized where it is applied, and 
the SFR dominates the rest of the flow. For example, 
consider annular Couette flow. As can be seen in Fig- 
ure ^1 the slip-lines only coincide with the shear planes 
(which in this case are concentric circles) along the inner 
wall of the cell. If the inner wall were given an increase 
Ar in applied shear stress, torque balancing requires the 
shear stress along any concentric circle within the mate- 
rial to receive a boost of Ar • Vw/r. Suppose the shear 
cell has inner wall radius 40(i and the boost in wall shear 
is significant, say Ar = r^^/lO. Solving for t — jicr along 
the shear planes in this situation gives a very different 
result than in the previous case — here r will only exceed 
ficr along the shear planes that are less than lAd off the 
inner wall. So, regardless of whether the density is or is 
not uniform, Bagnold scaling would at best only apply 
in an almost negligibly thin region near the wall. If the 
wall friction were less than fully rough, this region would 
further decrease. 



D. Redistribution of excess shear stress 

A more rigoruous physical justification of our conjec- 
ture can be made utilizing limit-state stresses and observ- 
ing the effect of pushing the system above yield. Bagnold 
rheology is a statement connecting the shear stress excess 
(i.e. amount by which r exceeds /i<j) along a shear plane 
to the rate of simple shear along the plane. We must em- 
phasize that shear planes and slip-lines are not equivalent 



E. A simple composite theory 

The preceding discussion indicates that, in general, one 
can use the admissibility status of the system to choose 
whether or not the flow should obey the SFR or Bag- 
nold rheology, or perhaps some combination of the two. 
Indeed, it seems reasonable that when slip-lines are ap- 
proaching admissibility (e.g. plate-dragging with high 
qo) or when an admissible system is only slightly pushed 
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above yield (e.g. inclined plane flow near static repose) 
we must account for contributions from both effects si- 
multaneously. It is beyond the scope of this paper to 
postulate the precise microscopic dynamics (and derive 
corresponding continuum equations) for this regime, but 
we can at least give a sense of how the more general the- 
ory might look. 

In general, we envision a smooth transition from rate- 
independent SFR dynamics to rate-dependent Bagnold 
dynamics controlled by the distribution of shear stress 
excess. This implies the coexistence of (at least) two 
different microscopic mechanisms: SFR and admissible 
shear. The SFR contribution would derive from the 
usual spot-based, quasi-static stochastic dynamics; the 
Bagnold contribution would come from a rate-dependent 
shearing motion along the appropriate slip-line family 
whenever there is a small excess shear stress (beyond the 
limit state) applied on a boundary which causes shear 
stress excess along the shear planes within. 

The two mechanisms should have different statistical 
signatures. For shear deformation along admissible slip- 
lines, we would expect anisotropic velocity correlations. 
In the direction perpendicular to the shear plane, the 
correlation length should be somewhat shorter than the 
typical spot size, since slip-line admissibility allows flow 
to occur with less drastic local rearrangements, farther 
from jamming. In the directions parallel to the shear 
plane, however the correlation length could be longer, 
since material slabs sliding along shear planes may de- 
velop more rigid, planar regions. It would certainly be 
interesting to study velocity correlations in heap flows at 
different inclinations and plate-dragging experiments un- 
der different loads to shed more light on the microscopic 
mechanisms involved in the SFR to Bagnold transition. 

For the remainder of this section, we make a first at- 
tempt at a composite model, simply a linear superposi- 
tion of SFR and Bagnold velocity fields: 

u = ausp^+/3u3^g. (51) 

which could have its microscopic basis in a random com- 
petition between the two mechanisms, when slip-lines are 
near admissibility. Here, u^^^ is an SFR solution for the 
flow, using the limit-state stress field everywhere, and 
Up is obtained from the excess shear stress on a bound- 

Bag 

ary by integrating the Bagnold strain rate 7 = — fiG 
over those shear planes for which r — /j.a has been boosted 
above zero. (Note that we ignore the condition of uni- 
form density for Bagnold rheology since we conjecture 
that uniform density is actually a geometric consequence 
of slip-line admissibility and will arise naturally whenever 
Bagnold rheology dominates the flow.) 

A reasonable first approximation is that the SFR 
and Bagnold solutions individually fulfill the necessary 
boundary conditions for the velocity profile since, under 
the right circumstances, either can be made to dominate 
the other. The constant [3 is the Bagnold proportional- 
ity constant which may depend on the density of the flow 
among other parameters 85]. Since the SFR is a rate- 



independent flow model, u^^^ can always be multiplied 
by a positive constant (observe that equation (|26|) is ho- 
mogeneous in ps), and thus we allow the scalar multiple 
a. Given some determinable form for a is chosen such 
that u fits the velocity boundary conditions. This seems 
reasonable for moving walls (as in plate dragging), but 
not for free surfaces, whose boundary velocity should also 
be predicted by the theory. In such cases, where a is not 
clearly defined in this simple model, one could use other 
empirical relations, such as the Pouliquen Flow Rule for 
inclined plane flows [s^, to deduce the free boundary 
velocity, and thus a. 

F. Some applications of the composite theory 

Using our very simple composite theory, we will now re- 
visit a few geometries that were troublesome for the SFR 
alone. Extending the theory with a smooth transition 
to Bagnold scaling controlled by slip-line admissibility 
seems to resolve the experimental puzzles and capture 
the basic physics of granular shear flows. In the cases 
we consider below, we do not change the value of a as 
we increase the shear stress excess; this way the relative 
importance of Bagnold effects are easier to isolate. 

1. Planar shear cell 

In a zero-gravity planar shear cell, u^^^ = 0, but the 
Bagnold solution for any amount of shear stress excess 
is of the form u^^^ = ky^ and thus the composite solu- 
tion, regardless of the values a and is a homogeneous 
flow between the two rough plates. The lack of a "back- 
ground" SFR solution in this case may explain why Bag- 
nold rheology is almost exactly observed in simulations 
of this geometry over a wide range of strain rates |6^ . 

2. Rough inclined plane 

For flow down a rough inclined plane, the SFR solution 
is an exponential decay (^TJ, and the Bagnold solution is 
a 3/2 power law (|5n|) . When the material is only slightly 
above static repose, a shear stress excess along the shear 
planes will exist but will be very small; it goes as V 
for an incline above static repose 86] . As a result, 
u„ will be small in mamitude, and the SFR solution 

Bag O ' 

will show through as the "creeping flow" with exponen- 
tial decay. As AO increases, the increased shear stress 
excess will cause the Bagnold contribution to increase, 
and the flow will eventually morph into the 3/2 power 
law dependence of pure Bagnold scaling. In between, 
where both contributions are of similar magnitude, the 
superposition of the two flow fields gives a predicted pro- 
file that appears approximately linear, since the SFR and 
Bagnold solutions are of opposite concavity. Thus, the 
composite SFR-Bagnold formulation appears to be able 
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FIG. 22: Some predicted velocity profiles for flow down a 
rough inclined plane as a function of depth. Note that all 
three known flow behaviors; exp decay, linear, 3/2 power 
law; appear in proper relationship to the inclination (see Ref. 
173). (Bottom) Incline near static repose, fully SFR domi- 
nated (Ls = 4:d); (Middle) Increased inclination angle, Bag- 
nold to SFR ratio of 3:1; (Top) Further increase to inclination; 
fully Bagnold dominated. 

to explain the various flow regimes in inclined plane flow, 
which have been observed in experiments and simulations 
(see Figure !^ . 

Recent experimental work of Pouliquen ^S§\ seems to 
support this analysis; it is found that inclined plane flow 
occurring at lower inclination angles exhibits spatial ve- 
locity correlations near the typical spot size (as the SFR 
would imply), but as inclination increases, the correlation 
length appears to decrease, an effect we might attribute 
to an increased dominance of Bagnold scaling (a phe- 
nomenon not goverened by a correlation length) over the 
SFR. 




FIG. 23: Standard heap flows enable one to see both the SFR 
and Bagnold contributions spearately in one flow geometry, 
(a) The surface region is dominated by Bagnold scaling, (b) 
The creep-region beneath adheres to the SFR. 
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FIG. 24: Plate-dragging slip-lines approaching admissibility 
as qo increases. 



3. Rapid heap flows 

The composite theory also seems consistent with rapid 
heap flows. When the flow rate down the heap increases, 
the region near the surface resembles inclined plane flow 
in any one of its various flow regimes, whereas the region 
beneath the surface flow undergoes creep motion which 
decays close to exponentially 33, 78] (see Figure We 
can justify this in terms of slip-line admissibility: The 
slip-lines throughout the system (see Figure [?T|) have the 
same form in both regions. In the surface region, the slip- 
lines are admissible because there is nothing blocking the 
motion from being a simple shearing along the slip-lines. 
In the creep region, however, the gate (or the ground) 
prevents global shearing along the slip-lines and thus the 
slip-lines are inadmissible and the SFR dominates. 

We can equivalently explain heap flow in terms of shear 
stress excess. The excess incurred by increasing the heap 
angle will distribute itself differently in the two regions. 



In the surface region, the excess can only be absorbed 
along the shear planes by inducing a strong Bagnold de- 
pendence. However, the gate at the bottom of the heap 
will support any shear stress excess on the creep region. 
(Note that the slip-lines in the creep region all hit the 
gate, or the ground.) Thus the full flow will be the sum 
of an exponentially decaying SFR solution superposed 
with a significant Bagnold-type solution which starts at 
the surface and cuts off at the interface with the creep 
zone. 



4- Plate dragging under a heavy load 

We will now explain the comment made at the end 
of the plate-dragging section. The slip-lines of a plate- 
dragging geometry can be pushed drastically close to full 
admissibility by simply increasing the relative loading of 
the top plate, go, above a certain non-excessive amount 
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(see Figure To see the effects of approaching admis- 
sibility more carefully, say we take a limit-state plate- 
dragging setup and pull the plate slightly harder, induc- 
ing a super-yield shear stress boost of Ar under the plate. 
Bagnold effects should appear wherever, as a consequence 
of stress balancing, a shear stress excess results along a 
shear plane. The shear planes are horizontal lines in this 
case, and at limit-state, the stresses along any horizontal 
obey 

r - /icr = -fifgV -j^—y- (52) 

When the extra shear stress is applied, a shear stress 
excess of Ar — /j^Pov/qo will form for < < -^Qo- 
Accordingly the Bagnold flow contribution will have a 
shear zone whose depth extends into the granular bed 
as an increasing linear function of qo for a fixed shear 
stress boost Ar and fixed downward plate-pressure (i.e. 
we decrease the material weight density to increase qo). 
Integrating the Bagnold shear rate gives that the relative 
size of the Bagnold contribution should also increase with 
increasing qo. With qo large enough, therefore, the SFR 
contribution will be dwarfed by a Bagnold term with a 
larger shear band. As g'o ^ oo the slip-lines become 
completely admissible and the shear band width diverges 
as we would expect. 

The mismatch in shear band size between the data 
of Tsai and Gollub and the SFR could be due, among 
other possible reasons, to the fact that qo ^ 3S0d was 
large enough to make the Bagnold contributution size- 
able. This could also explain the large shear bands found 
in [t^ . A detailed comparison of experiments and simula- 
tions with different versions of a composite SFR/Bagnold 
theory would be an interested direction for future work. 



VII. CONCLUSION 
A. Highlights of the present work 

We have proposed a stochastic fiow rule (SFR) for 
granular materials, assuming limit-state stresses from 
Mohr-Coulomb plasticity (MCP). In the usual case where 
slip-lines are inadmissible (inconsistent with boundary 
conditions), we postulate that fiow occurs in response to 
diffusing "spots" of local fiuidization, which perform ran- 
dom walks along slip-lines, biased by stress imbalances. 
The spot-based SFR corrects many shortcomings of clas- 
sical MCP and allows some of the first reasonable fiow 
profiles to be derived from limit-state stresses, which en- 
gineers have used for centuries to described the statics of 
granular materials. 

Our theory notably differs from all prior continuum 
theories (cited in the introduction) in that it is derived 
systematically from a microscopic statistical model p^ . 
The Spot Model is already known to produce realistic 
fiowing random packings 28], and what we have done is 



to provide a general mechanical theory of spot dynam- 
ics. Evidence for spots has been consistently found in 
spatial velocity correlations in simulations and ex- 
periments 57] (Fig. [TT|) on silo drainage. 

Beyond its fundamental physical appeal, the SFR 
seems to have unprecedented versatility in describing dif- 
ferent granular fiows. It has only two parameters, the 
friction angle (j) and correlation length (spot size) L^, 
which are not fitted; they are considered properties of 
the material which can be measured independently from 
fiow profiles. For monodisperse frictional spheres, the 
SFR can predict a variety of different fiows using the 
same spatial velocity correlation length, Lg ^ 3 — bd, 
measured in experiments and simulations. This is per- 
haps the most compelling evidence in favor of the spot 
mechanism which underlies the SFR. 

We have shown that the SFR can describe a rather di- 
verse set of experimental data on granular fiows. Some 
fiows are driven by body forces (silo and heap fiows); 
others have body forces, but are driven by applied shear 
(plate-dragging); still others are driven by applied shear 
without body forces (annular shear fiow). Some ge- 
ometries have straight boundaries (silos, heaps, plate- 
dragging), and yet the theory works equally well for 
highly curved boundaries (annular shear flow). Some of 
the flows exhibit shear localization (annular shear flow, 
plate-dragging, heap flow), and yet the theory correctly 
predicts wide shear zones in silo flow. It is notewor- 
thy that the same, simple model, correctly predicts and 
places shear bands in geometries where they arise for very 
different reasons — gravity causes the shear band in plate- 
dragging, and yet the geometry (through the Vp term in 
the drift) causes the shear band in annular Couette flow. 
We are not aware of any other theory (including classical 
MCP) which can quantitatively describe more than one 
of these flows, let alone without empirically fltting the 
velocity proflles. 

B. Comparison with partial fiuidization 

It is interesting to compare our approach to the contin- 
uum theory of partial fluidization of Aranson and Tsim- 
ring 31^ ^]. Although it lacks any microscopic basis, 
their theory also introduces a diffusing scalar fleld to con- 
trol the dynamics, as opposed to a classical stress/strain- 
rate relation. The analog of our spot density is the "order 
parameter" p, which measures the degree of "fluidiza- 
tion" of the continuum by mixing two different types of 
stresses, corresponding to distinct "liquid" {p = 0) and 
"solid" (p = 1) phases. Given the stress tensor for the 
material in a static solid state, a^^-, the stresses in a flow- 
ing granular material are modeled by adding some degree 
of viscous stresses, as in a Newtonian liquid: 

= (P + (1 - p)^iJ)^^J + r]E^J (53) 

where 77 is the viscosity. The order parameter controlling 
the balance of these two stress tensors is postulated to 
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obey a reaction-diffusion equation, 

{At)^ = l'V'p + p{l-p){p-S) (54) 

for collision time At, grain length scale and a function 
S of the stress state, which is greater than 1 where the 
material is above the static yield criterion, less than 
where below the dynamic yield criterion, and between 
and 1 otherwise. One benefit of this model is that it can 
be used for unsteady flows. In principle, the SFR may 
also describe time-dependence through the spot Fokker- 
Planck equation ((2^ . but we have only developed and 
tested the theory so far for steady flows, starting from 

For the sake of comparison, consider a steady flow 
modeled by partial fluidization and the SFR. The dif- 
ference is that the spot equation couples diffusion 
to a drift depending on frictional yielding, whereas the 
order parameter equation balances diffusion with a 
nonlinear source term, ressembling a chemical reaction 
rate, which indirectly mimicks the effect of a Coulomb 
yield criterion. Interestingly, if the SFR could be ex- 
tended to an elasto-plastic model without making the 
incipient failure assumption (see below), a similar non- 
linear source term may have to be added to the spot 
equation to account for the need to destroy spots when 
they enter zones below yield. It is also notable that our 
argument for why a spot drifts, i.e. a localized stick-slip 
type of shear stress decrease along the spot boundary, 
is reminiscent of equation wherein the shear stress 
goes down in the presence of fluidization. In this sense, 
a higher spot concentration in our model is similar to a 
higher degree of partial fluidization. 

The similarities, however, seem to end there. One dif- 
ficulty with the partial fluidization approach is that it 
cannot easily describe rate-independent effects since the 
motion stems from a viscous form in the stress tensor. 
Also in sharp contrast to our approach based on plas- 
ticity, partial fluidization does not provide a clear the- 
ory of the static solid stresses in the limit of no flow, 
opting instead to deal with environments for which the 
open components of this tensor are not needed (simple 
shear flows). This could perhaps be modified, but if the 
theory is to be fully general, equation (|53|) would need 
a frame-independent form where foreknowledge of the 
shear planes is not necessary to properly apply fluidiza- 
tion. These considerations as well as selecting boundary 
conditions on the order parameter, seem to be the pri- 
mary limitations in testing partial fluidization in more 
general situations, such as those considered here with the 
SFR. 



C. Future directions 

In spite of some successes, we still do not have a com- 
plete theory of dense granular flow. There are at least 
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FIG. 25: The difference between 2D and 3D continuum cells 
which could be used to construct the SFR. 

three basic limitations of the SFR, about which we can 
only offer some preliminary ideas to guide future work. 



1. Slip-line inadmissibility 

Although most slip-line fields are inadmissible, the 
SFR breaks down as slip-lines approach admissibility. 
We have already begun to extend the model into this 
regime by conjecturing that slip-line admissibility is as- 
sociated with Bagnold rheology, as excess shear stress 
(above the limit-state) drives a local shear rate along the 
most admissible slip-lines. We have shown that a sim- 
ple linear superposition of Bagnold and SFR flow fields 
with appropriate boundary conditions can describe a va- 
riety of composite flows, exhibiting both Bagnold and 
SFR behavior in different limits or segregated into dif- 
ferent regions. These include planar zero-gravity shear, 
various inclined-plane and heap flows, and plate-dragging 
at large loading. However, more work is needed to de- 
velop and test a composite SFR/Bagnold theory, both at 
the continuum level and in terms of the two microscopic 
mechanisms. 



2. 2D symmetry 

Through MCP, the SFR is currently used only in 
quasi-2D geometries. In efforts to extend the theory 
to 3D, a good test case would be the split-bottom Cou- 
ette cell, which displays a wide, diffusive shear band [H^, 
reminscent of a draining silo. The 2D limitation may not 
be so difficult to overcome, although any plasticity theory 
is more complicated in three dimensions, than in two. As 
usual, constructing a 3D limit-state stress field requires 
an additional hypothesis to close the stress equations. In 
3D, a general material point at incipient failure with dis- 
tinct principal stresses ai > (72 > os is intersected by a 
pair of slip planes angled 2e apart. We cannot therefore 
encase a 3D cell of material within slip planes as we are 
able to do to a 2D cell with slip-lines. However, the prin- 
cipal plane on which (72 acts, the intermediate principal 
plane, can be used along with the slip-planes to encase a 
3D material cell. This is legitimate because, if such a cell 
underwent slip-plane fluidization, the net material force 
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would be guaranteed to point parallel to the intermediate 
principal plane; since the intermediate principal plane of- 
fers no shear resistance, the material can slide along this 
plane, while simultaneously sliding along a slip-plane. 

To apply the SFR then, the drift vector should still 
be calculated from equations (|28j) and ((30j) . but all vec- 
tors must be projected first into the ai as-plane since 
the (72 direction is not involved in slip-plane fluidization. 
The shape of a spot and its diffusivity would likely be 
anisotropic, with different values in the intermediate di- 
rection, since the main source of diffusion is slip-plane 
fluidization. 

If ever the intermediate principal stress equals either 
the major or minor principal stress, as in the Har Von 
Karman hypothesis, incipient failure is upheld on a cone 
instead of intersecting slip-planes. When this degenerate 
case occurs, the material cell can be encased server al dif- 
ferent ways depending on the surrounding stress states. 
This must be determined before we can rigorously define 
how to apply fluidization and the SFR. 



3. Incipient yield everywhere 

The SFR assumes stresses near a limit state. While 
incipient yield is believed to be a good hypothesis in 
many situations, even during dense flow (which should 
be checked further in DEM simulations), it clearly breaks 
down in some cases, at least in certain regions. This is 
perhaps the most difficult limitation to overcome, since 
the limit-state assumption is needed to fully determine 
the stress tensor. Without it, the material effectively en- 
ters a different state most likely governed by some non- 
linear elastic stress strain law which is far more difficult 
to apply. 

We have already argued that such a transition away 
from incipient yield must exist in some granular flows, 
due to the strong tendency of granular materials to com- 
pactify into a rigid solid state when shaken (e.g. by 
nearby flowing regions) but not sufficiently sheared. A 
good example is a tall narrow silo with smooth side walls, 
where the SFR holds near the orifice, but breaks down 
in the upper region, ressembling a vertical chute. The 
broad shear band localizes on the side walls, as a rigid 
central plug develops, which likely falls below incipient 
yield. 

A more robust elasto-plastic theory for the stress state 
would relax our limit-state constraints and allow for ma- 
terial to fall below the yield criterion where it is described 
by elasticity. The SFR could then be applied only where 
the material is at yield and everywhere else the material 
does not deform plastically. Elasto-plasticity theory op- 
erates just as well in 3D as in 2D which is a key benefit 
over limit-state plasticity. However, our model as we have 
already presented it is far simpler than elasto-plasticity 
and yet still manages accurate results when applied to 
limit-state stress fields. 



As the SFR matures as a theory of granular flow, it 
would also be interesting to apply it to other amorphous 
materials, such as metallic glasses, and to develop new 
simulation methods. The basic idea is very general and 
applies to any material with a yield criterion. It has al- 
ready been suggested that the Spot Model could have 
relevance for glassy relaxation 27], and the SFR pro- 
vides a general means to drive spot dynamics, based on 
solid mechanical principles. The Spot Model also pro- 
vides a multiscale algorithm for random-packing dynam- 
ics, which works well for silo drainage [2^, so the SFR 
could enable a general framework for multiscale model- 
ing of amorphous materials. The idea would be to cycle 
between continuum stress calculations, meso-scale spot 
random walks, and microscopic particle dynamics. 
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APPENDIX A: CRITICAL STATE SOIL 
MECHANICS 

A common precept in plasticity is the notion of nor- 
mality or associatedness. Flows that obey normality have 
a flow rule defined in terms of the yield function Y as fol- 
lows: 



E = A— - 
dT 



(Al) 



where A is a positive multiplier. For a 3D flow, this means 
that if the yield function were plotted in 6-space as a func- 
tion of all 6 independent entries in the 3D stress tensor, 
the strain-rate matrix would be a 'vector' pointing nor- 
mal to the yield surface oriented toward greater values of 
Y. 

One of the first gripes about the use of friction-based 
yield criteria in describing granular materials is that the 
principle of normality gives a flow rule that predicts 
unstoppable dilatancy. Consider a rough extension of 
the Coulomb yield criterion into 3D, F = /i(trT)/3 
I To I /a/2, which displays the basic property that yield 
occurs when a certain multiple of the pressure equals the 
shear stress. Its associated flow rule is 



E = A 



72 1 To I 
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The trace of this strain-rate tensor is A/i, implying that 
material undergoing plastic flow will never stop dilating. 

Roscoe and co-workers present a different view- 
point on the issue. In what became known as Critical 
State Soil Mechanics, explained in detail in jsH, they 
argue that normality does hold, but that in fact the 
Coulomb yield criterion is not technically the correct 
yield function. 

Backed by results from triaxial stress experiments on 
soil samples, Critical State Theory claims that soils have 
a yield function that depends on the soil consolidation 
as measured by the local density p. The yield curve for 
material at a particular density is defined in terms of two 
stress tensor invariants: the pressure p = — ^trT and the 
equivalent shear stress q = |To|/a/2. Plotted in these 
variables, the principle of normality is equivalent to the 
statement that the strain-rate vector (e,7) is normal to 
the yield curve and pointing outward, where e = — trE 
is a volumetric strain-rate which determines changes in 
density, and 7 = ^|Eo| is a shear strain-rate propor- 
tional to the total shear deformation (volume-conserving 
part of the deformation). Figure [211 displays the theory's 
picture of the yield function and how it changes after 
material deformation. Any stress state underneath the 
yield curve corrresponds to rigid material. Under nor- 
mality, material at point 1 in the initial state will undergo 
a deformation according to the vector (ei,7i). Since ci 



is negative, the material will dilate and settle down at 
point 1 on the right on a new yield curve corresponding 
to < p. The material at stress state 2 will likewise 
undergo compaction and arrive at point 2 on the yield 
curve corresponding to > p. 

The critical state line is defined as the locus of points 
for which normality predicts no volumetric changes dur- 
ing deformation — note that wherever a yield curve inter- 
sects the critical state line, the curve becomes parallel to 
the p axis and thus the corresponding strain-rate vector 
has no volumetric component. The theory reasons that 
the critical state line is indeed a straight line of the form 

q = Mp. 

As flow developes, the stress states throughout the ma- 
terial will continually move toward the critical state line 
and once local volume changes finally stop, all flowing 
material stress states should lie on the critical state line. 
Thus in a steady flow, the critical state line might falsely 
appear to be the yield function when in fact it is only a 
locus of states from a family of yield functions. So, it is 
argued then that the reason normality previously failed 
to describe granular materials was because it was applied 
mistakenly to the critical state line and not to the true 
family of yield functions. 



R. M. Nedderman, Statics and Kinematics of Granular 
Materials (Nova Science, 1991). 

V. V. Sokolovskii, Statics of Granular Materials (Perga- 
mon/Oxford, 1965). 

R. Hill, The Mathematical Theory of Plasticity (Oxford 
at the Clarendon Press, 1950). 

H. M. Jaeger, S. R. Nagel, and R. P. Behringer, Rev. 
Mod. Phys. 68, 1259 (1996). 

P. G. de Gennes, Rev. Mod. Phys. 71, S374 (1999). 

L. P. Kadanoff, Rev. Mod. Phys. 71, 435 (1999). 

S. F. Edwards and D. V. Grinev, Advances in Complex 

Systems 4, 451 (2001), (also reprinted in Ref. |3). 

T. Halsey and A. Mehta, eds., Challenges in Granular 

Physics (World Scientific, 2002). 

I. S. Aranson and L. S. Tsimring, Rev. Mod. Phys. 78, 
641 (2006). 

C. S. Campbell, Ann. Rev. Fluid Mech. 22, 57 (1990). 
N. Balmforth and A. Provenzale, Geomorph. Fluid Mech. 
(Springer, 2001). 

J.-P. Bouchaud, M. E. Gates, J. R. Prakash, and S. F. 

Edwards, J. Phys. I (France) 4, 1383 (1994). 

J.-P. Bouchaud, M. E. Gates, J. R. Prakash, and S. F. 

Edwards, Phys. Rev. Lett. 74, 1982 (1995). 

T. Boutreux, E. Raphael, and P.-G. de Gennes, Phys. 

Rev. E 58, 4692 (1998). 

T. Boutreux, H. A. Makse, and P.-G. de Gennes, Euro. 
Phys. J. B 9, 105 (1999). 

R. A. Bagnold, Proc. Roy. Soc. London Ser. A 225 
(1954). 



[17] D. Ertas and T. C. Halsey, Europhys. Lett. 60, 931 
(2002). 

[18] S. B. Savage, J. Fluid Mech. 377, 1 (1998). 

[19] W. Losert, L. Bocquet, T. C. Lubensky, and J. P. Gollub, 

Phys. Rev. Lett. 85, 1428 (2000). 
[20] L. Bocquet, W. Losert, D. Schalk, T. C. Lubensky, and 

J. P. Gollub, Phys. Rev. E 65, 011307 (2002). 
[21] P. Mills, D. Loggia, and M. Tixier, Europhys. Lett. 45, 

733 (1999). 

[22] O. Pouliquen and R. Gutfraind, Phys. Rev. E 53, 552 
(1996). 

[23] O. Pouliquen, Y. Forterre, and S. L. Dizes, Advances 
in Complex Systems 4, 441 (2001), (also reprinted in 
Ref. 8]). 

[24] J. Litwiniszyn, Rheol. Acta 2/3, 146 (1958). 
[25] J. Mullins, J. Appl. Phys. 43, 665 (1972). 
[26] R. M. Nedderman and U. Tiiziin, Powder Technology 22, 
243 (1979). 

[27] M. Z. Bazant, Mechanics of Materials 38, 717 (2006). 
[28] C. H. Rycroft, M. Z. Bazant, G. S. Grest, and J. W. 

Landry, Physical Review E 73, 051306 (2006). 
[29] A. Lemaitre, Phys. Rev. Lett. 89, 195503 (2002). 
[30] A. Lemaitre, Phys. Rev. Lett. 89, 064303 (2002). 
[31] L S. Aranson and L. S. Tsimring, Phys. Rev. E 64, 

020301 (2001). 

[32] L S. Aranson and L. S. Tsimring, Phys. Rev. E 65, 
061303 (2002). 

[33] G. D. R. Midi, Euro. Phys. Journ. E. 14, 341 (2004). 
[34] M. Ostoja-Starzewski, Int. J. Plasticity 21, 1119 (2005). 



29 



> 

q 


/ Critical State Line 




> 










// Rigid stres^^^^^ 
//' states 














P 


P 



FIG. 26: (left) Critical State Theory's spade-shaped yield function for material at some density p^; (right) Deformations with 
non-zero volumetric strain will cause the material to settle down on a new yield function. 
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